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FOREWORD 


This  Technical  Documentary  Report  has  been  prepared  in 
four  volumes,  as  follows: 

Contractor 1 s 

Volume  _ Title _  Publication  Number 

I  Program  Development  U-  3005 

II  Operating  Instructions  U-  3006 

III  Programmer's  Manual  U-  3007 

IV  Operations  Summary  (U)  S-  2990 


Publication  of  this  technical  documentary  report  does 
not  constitute  Air  Force  approval  of  its  findings  or 
conclusions.  It  is  published  only  for  the  exchange 
and  stimulation  of  ideas. 


SPIRAL  DECAY  AND  SENSOR  CALIBRATION 
DIFFERENTIAL  CORRECTION  PROGRAMS 


ABSTRACT  (U) 


This  document  presents  the  theory,  operation,  and  coding  details 
of  and  experience  with  a  computer  program  for  the  accurate  prediction  of 
the  reentry  corridor  of  an  earth  satellite  undergoing  orbital  decay  due  to 
atmospheric  friction.  This  program  is  also  useful  in  the  precision  pre¬ 
diction  of  satellite  positions  for  other  purposes  such  as  sensor  calibration 
and  weapon  systems. 
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SECTION  I 

INTRODUCTION 


This  report  documents  Aeronutronic 1 s  approach  to  the  formulation  of  a 
computer  program  for  high-precision  orbit  determination  within  the  SPACETRACK 
computational  environment.  Although  originally  intended  for  terminal  decay 
corridor  determination  for  satellites  executing  spiral  decay,  this  program  has 
evolved,  within  the  limitations  imposed  by  the  Spiral  Decay  Technique  Development 
Contract,  as  a  highly  effective  tool  for  a  majority  of  the  high-precision  require¬ 
ments  of  the  SPACETRACK  center,  including  spiral  decay,  weapon  command  and  control, 
and  sensor  calibration.  Through  an  extended  period  of  joint  contractor  and 
SPACETRACK  analyst  application  of  this  class  of  programs  to  the  rapidly-changing 
spectrum  of  operational  requirements  on  the  SPACETRACK  center,  a  high  degree  of 
analyst  proficiency  and  operational  utility  have  also  been  developed. 

Specifications  for  the  Spiral  Decay  class  of  programs  are  conveyed  in 
Table  I.  Two  programs  are  involved: 

(1)  Spiral  Decay  (SPIRDEC)  is  an  orbit  determination  and 

prediction  program,  with  weighted  differential  correction 
capability  for  osculating  elements  and  ballistic  coeffi¬ 
cient.  Current  limitations  are  to  non-equatorial 
satellites  with  eccentricities  less  than  0.9*  (no 
ballistic  coefficient  correction)  and  0.1  (with  ballistic 
coefficient  correction) . 


*  For  the  higher  eccentricities,  lunar  and  solar  perturbations  become  increasingly 
important.  These  are  not  included  in  the  current  version  of  the  program. 
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TABLE  I.  SPIRAL  DECAY  PROGRAM  SPECIFICATIONS 


Operating  System 

SPS  B-3  /  Philco  2000 

Schedule  Tape  Mode 

Ephemeris  Integration 

Variation  of  Parameters 

Cowell  Entry  Option 

Numerical  Integration 

Sixth-order  Adams-Bashf or th  with 
Runge-Kutta  Start 

Automatic  Interval  Control 

Initialization 

SPACETRACK  Mean  Elements 

Osculating  (M,  N)  Elements 

Geocentric  Position  and  Velocity 

Geopotential 

Zonal  Harmonics  Inclusive  of  Fifth- 
order 

Non-zonal  Harmonics  Inclusive  of 
Fourth-order 

Atmosphere 

Jacchia-Nicolet  Dynamic  Model  Above 

120  Kilometers 

COESA  1962  Static  Model  Below  120 
Kilometers 

Differential  Correction 
(Weighted) 

Six  Osculating  Parameters  (M,  N 
Formulation) 

Ballistic  Coefficient 

Sensor  Location  Timing  and  Observation 
Biases (Using  CALIB  Program) 

Sensor  Input  Data 

Internally  Compiled  Table  for  Biases  and 
Weights,  with  Optional  Override 

Observation  Capacity 

Up  to  984  Observations;  Sorting  Internal 
to  Program 

Prediction  Capability 

Geographic  and/or  Geocentric  Cartesian 
Coordinates 

Acquisition  Capability 

Binary  Ephermis  Tape  for  XYZLA  Program 

Restart  Capability 

Punches  Parameter  Cards  for  Both  Spiral 
Decay  and  Calibration  Programs 
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(2)  Sensor  Calibration  (CALIB)  is  a  differential  correction 

program  for  determining  biases  in  sensor  location,  timing, 
and  observed  quantities,  utilizing  a  reference  orbit 
determined  from  an  independent  sensor  network  (normally 
either  PRELORT  or  Baker-Nunn  data),  using  the  SPIRDEC 
program. 

The  SPIRDEC  program  has  been  compiled  in  both  SPS  B-2  and  SPS  B-3  versions  to 
accommodate  the  SPACETRACK  system  change  to  be  installed  during  January  1965, 
Although  the  operating  instructions  conveyed  with  this  report  reflect  a  B-3 
capability  for  the  CALIB  program  as  well,  delivery  of  a  B-3  version  coincident 
with  this  report  is  not  planned.  In  the  interim,  the  SPS  A-l  version  of  the 
SPIRDEC  and  CALIB  programs,  documented  in  the  Milestone  II  report, x  shall  be 
available . 


In  the  development  of  the  Spiral  Decay  class  of  programs,  particular 
attention  has  been  paid  to  those  factors  considered  limiting  in  a  useful  spiral 
decay  and  weapon  command  and  control  capability.  Important  considerations  which 
have  influenced  the  program  design  include: 

(1)  Environmental  factors  which  influence  the  trajectory, 
primarily  the  atmospheric  density  and  gravitational 
perturbations,  must  be  carefully  considered.  Such 
considerations  lead  to  the  inclusion  of  a  dynamic  model 
atmosphere  and  longitude-dependent  (tesseral)  harmonics 
in  the  geopotential,  inclusive  of  fourth-order  .'oc 

(2)  The  rational  application  of  statistical  weights  to  ADC 
sensor  data  cannot  be  made  without  first  providing  for 
the  large  known  biases  in  these  sensors.  The  Spiral 
Decay  program  will  accept  and  process  both  weight  and 
bias  data  for  sensors  and,  through  the  interfacing 
Sensor  Calibration  program,  will  determine  location, 
timing  and  observation  biases  for  selected  sensors, 
utilizing  either  ADC  or  independent  (PRELORT,  Baker- 
Nunn)  sensors  as  a  reference  network. 


*  nSpiral  Decay  and  Sensor  Calibration  Differential  Correction  Programs  - 
Milestone  II  Draft”,  31  March  1964,  Aeronutronic  Publication  U-2559. 

**  The  fifth  zonal  harmonic  is  also  included.  A  subroutine  is  now  in  checkout 
which  extends  the  potential  to  the  ninth-order  zonal  and  sixth-order  non- 
zonal  harmonics.  A  change  note  to  this  documentation  will  be  issued  upon 
completion. 
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(3)  Since  this  class  of  programs  will  normally  provide  nearly 
real-time  support  to  SPACETRACK  spiral-decay  and  weapon- 
oriented  requirements,  every  effort  has  been  expended  to 
provide  a  precision  orbit  capability  of  the  highest 
computational  efficiency.  This  is  reflected  in  the  choice 
of  the  variation-of-parameters  method  for  ephemeris  calcu¬ 
lation,  and  in  the  derivation  and  application  of  analytical 
differential  correction  equations."' 

To  accomplish  these  ultimate  objectives  while  providing  a  meaningful 
interim  capability,  a  series  of  programming  milestones  was  established.  These 
milestones,  with  technical  and  programming  targets,  are  listed  below: 


Milestone 


Technical  and  Programming  Target 


I  (31  December  1963) 


II  (31  March  1964) 

III  (30  June  1964) 

IV  (30  September  1964) 

V  (31  December  1964) 


Demonstrate  the  application  of  the 
variation-of-parameters  approach , 
with  analytic  correction  equations, 
to  a  decaying  satellite  (object  292) 

Demonstrate,  in  the  live  SPACETRACK 
computation  and  data  environment,  an 
interim  capability  for  Spiral  Decay 
calculations  (objects  632,  707) 

Deliver  an  interim  operational  (SPS  B-2) 
capability  for  Spiral  Decay  calculations, 
Demonstrate  an  interim  sensor  calibra¬ 
tion  capability  (objects  759,  811) 

Demonstrate  a  weapon  support  capability 
in  the  SPACETRACK  environment. 

Deliver  "final11  operational  (SPS  B-3) 
capability  for  Spiral  Decay  and  weapon 


*  This  approach  has  been  demonstrated,  in  the  SPACETRACK  environment,  to  exceed 
the  computational  speed  of  the  competitive  COWELL  approach  by  a  factor  of  four 
or  more,  with  better  accuracy. 
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Highlights  in  the  development  of  these  capabilities  include  the  scheduled  visual 
confirmation  of  two  satellite  decays  during  the  Milestone  II  demonstration 
(See  Section  4)  and  the  subsequent  scheduling  of  four  additional  visual  confirma¬ 
tions  . 


This  report  documents  the  Milestone  V  program  status.  Subsequent  sections 
and  appendices  provide  a  detailed  technical  description  of  the  Spiral  Decay  and 
Sensor  Calibration  programs. 
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SECTION  2 


ORBIT  DETERMINATION  PROGRAM 


Except  in  the  somewhat  trivial  case  of  a  point  mass  moving  in  a  pure 
central  force-field,  the  so-called  ,fTwo  Body  Problem, 11  it  is  generally  impossible 
to  derive  orbital  elements  from  observations  by  direct  analytical  means.  The 
process  usually  employed  is  to  create  a  dynamical  model  of  the  motion  and  then 
to  adjust  the  constants  introduced  into  this  model  until  the  best  possible  fit 
to  the  observations  is  obtained  in  the  sense  of  least  squares.  The  problem  is 
very  similar  to  the  one  faced  in  classical  celestial  mechanics;  the  solutions 
discussed  here  draw  heavily  from  the  theory  and  practice  of  celestial  mechanics. 

The  principal  differences  between  modern  orbital  theory  of  artificial 
satellites  and  classical  celestial  mechanics  theories  arise  in  part  from  the  fact 
that  the  artificial  satellite  environment  includes  perturbative  forces  such  as 
aerodynamic  drag,  earth-bulge,  and  the  like,  and  in  part  from  the  availability  of 
high-speed  computing  equipment  and  electronic  instrumentation  providing  range  and 
range-rate  data,  as  well  as  angles. 

2.1  DIFFERENTIAL  CORRECTION  PROCEDURE 

Reduced  to  its  simplest  form,  the  differential  correction  technique 
involves  the  following  five  steps  (identified  in  Figure  1): 

Step  1.  Raw  observation  data  is  corrected  for  known  systematic 
errors  (biases)  as  determined  from  calibration  efforts. 
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Differential  Corrections  To: 


Step 


Step 

0 


Six  Orbit  Parameters 
Ballistic  Drag  Parameter 


Observation 

Weights 


£) 

Step 


FIGURE  1.  SPIRAL  DECAY  PROCRAM  FLOW  CHART 


Step  2.  Beginning  with  a  tentative  estimate  of  the  injection  conditions 
and  other  "constants"  entering  into  the  equations  of  motion,  an 
ephemeris  is  produced  by  integrating  these  equations.  At  each 
time  for  which  an  observation  exists,  this  ephemeris  is  trans¬ 
lated  into  a  "representation"  for  those  observations  available 
from  the  particular  sensor. 

Step  3.  In  the  absence  of  unknown  observation  errors,  and  for  exact 

knowledge  of  the  injection  conditions  and  environment  governing 
the  motion,  the  representations  and  corrected  observations 
would  be  identical.  In  practice,  however,  they  are  not  identical, 
and  their  differences,  or  "residuals,"  serve  as  a  basis  for 
differentially  correcting  the  initial  constants,  as  shown  in  the 
following  steps  4  and  5. 

Step  4.  Cause  and  effect  relations  between  injection  conditions  and 

subsequent  observations  are  highly  nonlinear.  Where  corrections 
to  these  injection  conditions  and  subsequent  changes  in  residuals 
are  small,  the  problem  is  readily  linearized.  For  each  observed 
quantity,  a  scalar  differential  expression  relating  the  known 
residual  to  the  unknown  changes  in  injection  conditions  and  other 
parameters  defining  the  motion  is  derived. 


Step  5.  The  system  of  linear  equations,  one  for  each  observed  quantity 
and  normally  heavily  overdetermined,  is  solved  in  the  sense  of 
weighted  least  squares  for  the  differential  corrections  to  the 
injection  parameters.  Some  insight  into  the  convergence  of  the 
procedure  is  available  by  examing  the  weighted  sum-of-squares 
of  the  residuals,  and  the  variance-covariance  matrix  of  the 
differential  corrections. 


Normally,  steps  2  through  5  will  be  repeated  until  the  weighted  sum-of-squares  of 
the  residuals  becomes  stationary.  During  this  iterative  process,  observations 
which  do  not  fit  the  orbit  may  be  rejected. 

Execution  of  the  foregoing  five  basic  operations  is  quite  involved  and 
requires  meticulous  attention  to  detail.  Such  items  as  the  atmospheric  model, 
geopotential,  and  choice  of  reference  systems  must  be  correctly  handled  to  ensure 
the  greatest  precision. 
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2.2  PROGRAMMING  SYSTEM 


To  facilitate  the  eventual  application  of  the  Spiral  Decay  Program  within 
the  ADC  Spacetrack  Center,  yet  permit  interim  experimentation  and  simulation  at 
496-L  (Bedford,  Massachusetts),  Aeronutronic  and  the  ADC  Spacetrack  Center,  the 
Semiautomatic  Programming  System  (SPS) ,  A-l  version,  was  utilized  in  the  first 
milestone.  The  successful  demonstration  of  the  Milestone  II  version  in  the  ADC 
environment,  and  the  desire  on  the  part  of  ADC  personnel  to  utilize  these  procedures, 
made  an  interruptable  (B-2)  version  of  this  program  desirable.  Thus,  the  Milestone 
II  program  was  prepared  in  two  versions: 

(1)  An  interruptable  (B-2)  version  of  the  Spiral  Decay  module, 
permitting  application  of  this  procedure  during  periods 
where  the  ADC  computer (s)  are  providing  backup  for  other 
ADC  functions. 

(2)  A  noninterruptable  (A-l)  version  of  the  Spiral  Decay  and 
Sensor  Calibration  modules,  in  which  the  executive  system 
has  been  modified  to  prevent  truncation  of  time.  This 
latter  feature  is  considered  essential  to  meaningful  cali¬ 
bration  work. 

The  operating  instructions  for  both  A-l  and  B-2  program  versions,  operating  in  the 
Schedule  Tape  Mode,  were  conveyed  in  the  Milestone  II  draft  of  the  present  report. v 

The  present  report  applies  to  programs  operating  under  the  control  of  the 
B-3  system,  which  is  interruptable  and  uses  new  formats  accommodating  five-digit 
satellite  numbers. 

2.3  EPHEMERIS  INTEGRATION 

The  ephemeris  integration  represents  a  major  portion  of  the  computational 
burden  for  spiral  decay  calculations.  Since  no  satisfactory  theory  is  yet  available 
to  perform  this  integration  by  general  perturbations  (development  of  the  perturba¬ 
tions  into  series  and  integrating  term-by-term)  for  orbits  highly  perturbed  by  drag, 
the  special  perturbations  technique  (numerical  step-by-step  integration)  has  been 
employed . 


*  Aeronutronic  Publication  U-2559,  31  March  1964 


9 


The  coordinates  are  referred  to  the  mean  equinox  of  a  fixed  epoch 
and  the  true  equator  of  date.  For  detail  refer  to  the  definitions  of 
I,  J,  K,  etc.,  in  Appendix  VII. 

Of  the  several  formulations  available  in  special  perturbations, 
the  variation-of-parameters  method  has  been  employed  for  two  reasons: 

(1)  The  technique  supresses  the  dominant  central 
gravitation  term  in  the  geopotential  by  inte¬ 
grating  parameters  which,  in  the  absence  of 
perturbative  effects,  are  constant. 

(2)  The  parameters  integrated  in  the  variation-of - 
parameters  method  may  be  related,  through  scalar 
differential  expressions,  to  observed  coordinates 

at  future  times,  thereby  permitting  a  fully  analyti¬ 
cal  development  of  the  differential  correction. 

These  factors  contribute  to  computational  efficiency  by  (1)  permitting 
large  integration  step  size  (or  fewer  derivative  evaluations  per  revolution) 
and  (2)  facilitating  communication  between  the  ephemeris  and  differential 
correction  modules  of  the  program. 

Of  the  several  available  formulations  for  variation-of -parameters , 
the  ci,  h  formulation*  has  been  used,  since  it  exhibits  no  zero-eccentricity 
singularities. 

During  the  last  revolution  of  a  satellite  executing  spiral  decay, 
the  perturbative  accelerations  become  too  large  to  be  efficiently  accommodated 
by  the  var iation-of-parameters  method.  A  transition  to  the  Cowell  formu¬ 
lation  of  the  equations  of  motion  during  this  revolution  is  provided. 

a .  Numerical  Integration  Method 


The  variation-of -parameters  method  leads  to  seven  first-order 
differential  equations , with  time  or  one  of  the  anomalies  as  the  independent 
variable. 


*  A  complete  derivation  of  the  equations  is  given  in  Appendix  VI.  The 
mathematical  steps  involved  in  both  ephemeris  integration  and  differential 
correction  are  conveyed  in  Appendix  I. 
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Prior  to  Milestone  IV,  expediency  dictated  program  compilation  with  the 
Runge-Kutta  integration  subroutine,  a  fourth-order  technique  requiring  four  (4) 
derivative  evaluations  per  step.  This  self-starting  technique  permits  great 
flexibility  in  specification  of  time  interval,  but  at  the  expense  of  running  time. 
Comparable  precision  and  stability,  however,  is  available  with  higher-order  pre¬ 
dictor-corrector  methods  requiring  fewer  derivative  evaluations  per  revolution. 

The  final  program  thus  is  compiled  with  such  an  alternative  numerical  integration 
subroutine,  Adams-Bashf or th .  The  Runge-Kutta  technique  will  continue  to  be  used 
for  starting  purposes. 

Representation  of  observations  for  residual  computation  demands  that 
satellite  coordinates  at  arbitrary  observation  times  be  calculated.  Initially, 
and  again  for  expediency,  the  program  numerically  integrates  to  each  observation 
time.  For  situations  where  a  large  number  of  observations  exist,  such  as  sensor 
calibration  work  utilizing  raw  data  derived  from  the  AF  Satellite  Control 
Facilities,  the  additional  integration  step  for  each  observation  becomes  in¬ 
efficient.  To  avoid  this  unnecessary  computational  burden,  a  fifth-order  divided 
difference  interpolation  procedure  has  been  included,  permitting  the  interpolation 
for  observation  times  in  an  ephemeris  table  of  arbitrary  time  interval  (as 
determined  by  the  integration  error  control) . 

b .  Drag  Perturbations 

A  number  of  atmospheric  models,  which  include  dynamic  terms  to  represent 
deviations  from  "standard11  models,  are  in  use.  Those  dynamic  effects  which  have 
been  identified  from  satellite  acceleration  studies  include  seasonal  and  longitude- 
dependent  terms,  as  well  as  those  correlated  with  solar  radiation  and  planetary 
magnetic  effects. 

At  altitudes  below  300  km.,  representing  the  region  of  interest  during 
the  last  few  days  of  spiral  decay,  the  variations  in  density  due  to  important 
dynamic  terms  are  given  in  Table  II. 


11 


TABLE  II.  VARIABILITY  OF  DYNAMIC  ATMOSPHERE  TERMS* ** 
(orders  of  magnitude) 


Altitude 

Diurnal 

Seasonal 

Solar  (F1q) 

Magnetic (A^) 

150  km 

0 

0 

0.14 

0 

200  km 

0 

0.07 

0.54 

0.05 

250  km 

0.06 

0.10 

0.86 

0.09 

300  km 

0.14 

0.12 

1.00 

0.14 

At  higher  altitudes,  these  effects  are  even  greater. 

These  considerations  have  led  to  the  incorporation  of  a  dynamic  model 
atmosphere  in  the  spiral  decay  program.  The  Jacchia-Nicolet  model™  has  been 
utilized;  this  model  includes  significant  geometric  and  solar  radiation  terms 
above  120  km.  The  formulation  of  the  Jacchia-Nicolet  atmosphere  subroutine  is 
given  in  Appendix  II. 

c .  Gravitational  Model 


In  the  selection  of  terms  to  be  included  in  the  gravitational  model  for 
the  earth,  careful  consideration  must  be  given  to  terms  which  either  (1)  lead  to 
secular  and/or  long  period  terms  which,  if  neglected,  would  be  aliased  into  other 
effects,  or  (2)  can  lead  to  short  period  radial  displacements  with  consequential 
drag  coupling. 

Based  upon  these  considerations,  zonal,  sectorial,  and  tesseral  terms 
have  been  included  in  the  gravitational  model.  For  Milestones  I  and  II,  the  J ^ 
term  was  implemented.  A  general  expression  for  the  earth  model  permitting  the 
inclusion  of  a  consistent  set  of  the  harmonics  through  n,m  =  4,4  is  now  used. 

The  formulation  of  this  earth  model  is  conveyed  in  Appendix  I.  The  values  of  the 
coefficients  of  the  non-zonal  harmonics,  which  are  still  undergoing  redefinition 
by  several  investigators,  are  specified  by  two  parameter  cards  in  the  program 
input . 


*  H.  K.  Paetzold,  ffModel  for  the  Variability  of  the  Terrestrial  Atmosphere  above 
150  km",  submitted  to  COSPAR  working  group  IV,  "Reference  Atmosphere",  1962 
January . 

**L.  G.  Jacchia,  "Temperature  Above  the  Thermopause" ,  Smithsonian  Inst. 
Astrophysical  Observatory  Special  Report  No.  150,  April  22,  1964. 
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2.4  DIFFERENTIAL  CORRECTION 


The  differential  correction  involves  the  formulation  of  scalar  differ¬ 
ential  expressions  relating  observation  residuals  to  corrections  to  the  parameters 
defining  the  representation;  these  parameters  normally  include  six  initial  condi¬ 
tions  or  elements  defining  the  orbit  and  a  seventh  parameter  related  to  the 
interaction  of  the  satellite  with  drag  environment.  Where  sufficient  observation 
redundancy  exists  to  permit  the  determination  of  additional  parameters,  sensor 
biases  (range,  angles,  timing)  or  geometric  factors  entering  into  the  representation 
(such  as  geocentric  coordinates  of  the  observer)  may  also  be  determined. 

The  scalar  differential  expressions  take  the  following  form: 

n  BO 

AO,  =  Z  _ L_  AP. 

j=l  *P. 

where  0.  are  observed  coordinates  and  the  p.  are  parameters  to  be  corrected.  Two 
conditions  must  be  met  to  solve  for  the  APj>  namely: 

(1)  There  must  be  at  least  as  many  observations  as  parameters 
(i>n)  to  solve  the  system  of  linear  equations,  and 

(2)  The  selected  parameters  must  be  capable  of  being  related 

to  observed  coordinates  by  linear  differential  expressions, 
without  singularities  at  critical  parameter  values. 

The  first  condition  is  readily  met  with  modern  electronic  instrumentation,  and  the 
overdetermined  system  of  equations  is  solved  in  the  sense  of  weighted  least  squares. 
This  approach  assumes  that  observation  errors  are  statistically  uncorrelated,  a 
reasonable  assumption  for  SPACETRACK  observations  if  the  biases  can  be  predetermined. 
Under  these  conditions,  the  variance-covariance  matrix  of  orbit  parameters  reveals 
the  quality  of  the  determination  and  may  be  printed  out  for  inspection. 


The  second  condition  above  implies  that  the  partial  derivatives  exist  for 
the  complete  range  of  parameter  values  encountered  in  the  differential  correction. 
Attempts  to  differentially  relate,  for  example,  the  orbit  eccentricity  and  perigee 
argument  to  observed  coordinates  through  linear  differential  expressions  must  fail 
at  zero  eccentricity  for  the  eccentricity  J.s,  by  definition,  non-negative  and  the 
perigee  argument  indeterminate.  A  number  of  parameter  sets  exist  which  can  be 
successfully  related  to  observed  coordinates  through  linear  differential  expressions 
at  this  critical  eccentricity*.  Some  of  these  sets,  such  as  the  M,  N  pattern 


*  See,  for  example,  Aeronutronic  Report  U-880. 
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utilized  in  the  Spiral  Decay  Program,  may  be  differentiated  to  provide  a  completely 
analytical  formulation  for  the  differential  correction,  thereby  leading  to  sub¬ 
stantial  economies  in  computational  time. 

With  the  exception  of  the  drag  parameter,  the  M,  N  form  for  the  differ¬ 
ential  correction  has  been  adequately  documented*.  Analytical  differential 
expressions  for  the  drag  parameter  have  also  been  derived,  and  this  derivation  is 
included  as  AppendixIV  to  this  report. 

In  addition  to  the  differential  correction  features  cited  above,  the 
Milestone  III  version  was  to  have  included  limited  differential  correction  in  the 
Cowell  entry  option.  It  has  been  demonstrated  (see  Section  6),  however,  that  the 
variation-of-parameters  method  will  accommodate  tracking  data  extending  into  the 
incendiary  region  of  decay.  Experimentation  has  also  established  the  desirability 
of  including  geometric  terms  in  the  drag  correction  equations  arising  from  earth 
flattening,  which  can  lead  to  density  variations  equivalent  to  a  scale  height. 

The  formulation  of  these  expressions  is  given  in  Appendix  IV. 


*In  addition  to  U-880,  sources  include  the  Proceedings  of  the  JPL  Seminar  on 
Tracking  Programs  and  Orbit  Determination,  22-26  February  1960,  and  Koelle, 
ffHandbook  of  As tronautical  Engineering'*  McGraw-Hill,  1961,  Section  8.112. 
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SECTION  3 


SENSOR  CALIBRATION  PROGRAM 


In  Appendix  V,  scalar  differential  expressions  have  been  developed  which 
relate  observation  residuals  —  in  range,  range-rate,  azimuth  and  elevation  —  to 
biases  in  sensor  location,  time  reference  and  observed  coordinates.  This  theory 
has  been  implemented  in  a  Sensor  Calibration  Program,  designated  CALIB,  operating 
in  the  Schedule  Tape  Mode,  A-l  system.  This  Section  presents  a  description  of  the 
procedure;  operating  instructions  are  given  in  Volume  II.  Results  of  operational 
applications  of  these  techniques  are  conveyed  in  the  Operations  Summary,  Volume  IV. 


The  calibration  procedure  involves  the  consecutive  application  of  two 
differential  correction  procedures  as  shown  in  Figure  2. 

The  first  step  involves  the  definition  of  a  highly  precise  reference 
orbit  derived  by  processing  tracking  data  from  sources  of  high  confidence  level 
and  redundancy.  The  most  appropriate  tracking  system  for  this  purpose  is  the 
network  of  PRELORT  radars  operated  by  the  6594th  Aerospace  Test  Wing  (Sunnyvale, 
California)  in  support  of  AFSC  satellite  programs.  These  radars  track  a  co¬ 
operative  S-band  beacon  carried  by  most  AFSC  satellites,  and  the  tracking  data 
rate  for  this  network  is  in  excess  of  1500  observations/day  even  for  relatively 
low  (100  n.m.)  satellites.  The  quality  of  these  data  and  of  the  earth  and 
atmospheric  model  used  to  define  the  satellite  motion  is  demonstrated  by  the 
ability  to  fit  these  data,  even  over  relatively  long  time  spans,  to  the  order  of 
100  meters.  To  facilitate  the  utilization  of  these  data  within  the  ABC  computa¬ 
tion  environment.  AFSC  has  contracted  for  preparation  of  a  160A  computer  program 
to  translate  PRELORT  data  tapes  into  SPADAT  BCD  card  format,  and  to  prepare  a 
standard  teletype  tape  for  transmission  into  SPACETRACK. 
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REFERENCE  DATA 
(PRELORT,  B-N,  etc.) 


CANDIDATE  DATA 
FOR  CALIBRATION 


FIGURE  2.  SENSOR  CALIBRATION  PROCEDURE 


Additional  data  sources  can  be  exploited  for  performing  this  calibration 
function,  either  to  complement  the  FRELORT  data  or  to  replace  it  where  the  satellite 
does  not  carry  an  active  S-band  beacon.  Baker -Nunn  camera  data,  for  example,  lacks 
the  redundancy  and  intrinsic  metric  desirable  for  calibration  exercises,  but  will 
strengthen  the  PRELORT  data,  where  greatest  emphasis  is  placed  upon  range  data. 

The  design  of  the  Spiral  Decay  differential  correction  program  provides 
for  the  definition  of  the  reference  orbit  to  the  desired  precision.  The  necessary 
interface  between  the  SPIRDEC  and  CALIB  programs,  as  shown  in  Figure  2,  is  a 
vector  defining  the  parameters  of  the  reference  orbit;  provision  has  been  made  in 
the  SPIRDEC  program  to  format  and  print  the  parameter  card  images  for  control  of 
CALIB  directly. 

The  second  step  in  the  Sensor  Calibration  procedure  also  involves  a 
differential  correction  procedure.  Beginning  with  the  reference  ephemeris  and  with 
observations  derived  from  the  sensor  undergoing  calibration,  a  weighted  least- 
squares  determination  is  made  of  those  sensor  properties  —  sensor  location,  time, 
and  observation  biases  —  which  lead  to  the  best  fit  of  the  tracking  data  to  the 
reference  ephemeris.  This  step  is  performed  by  the  CALIB  differential  correction 
program,  which  permits  the  selective  determination  of  any  combination  of  the  seven 
sensor  properties,  data  permitting.  In  addition  to  these  biases,  CALIB  also 
determines  RMS  values  for  the  individual  observation  types,  to  be  used  as  weights 
in  subsequent  applications  of  the  sensor  data. 

Although  these  procedures  could  have  been  assembled  as  a  single  program, 
the  following  considerations  influenced  the  decision  to  prepare  two  programs  with 
a  simple,  program-defined  interface: 

(1)  More  flexible  use  of  memory  is  permitted,  particularly  in 
subsequent  B-2  and  B-3  versions  of  these  programs,  where 
certain  memory  regions  are  denied. 

(2)  The  simultaneous  determination  of  biases  and  orbit  parameters, 
in  a  single  pass  involving  both  PRELORT  and  ADC  sensor  data, 
may  alias  the  reference  orbit  parameters  unless  the  freedom 

to  determine  all  biases  in  the  candidate  sensor  is  permitted. 

Few  such  data  situations  will  exist. 

It  is  important  to  note  that  the  earth  model  and  atmosphere  constants  must  be 
identical  in  the  two  programs,  to  permit  the  accurate  reconstruction  of  the 
reference  ephemeris  in  the  CALIB  program. 
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The  current  operational  version  of  the  calibration  procedure  operates  in 
the  Schedule  Tape  Mode,  A-l  version.  Restriction  to  the  A-l  operating  system  was 
necessitated  by  the  manner  in  which  time  is  processed  in  executive  modules  of  the 
B-2  system,  where  the  representation  of  time  since  1950.0  as  a  floating  point 
number  introduces  errors  in  tens  of  milliseconds.  In  the  A-l  master  tape  utilized 
at  ADC/496L  for  calibration  work,  these  modules  have  been  modified  to  represent 
time  in  fixed  point  and  double  precision,  thereby  avoiding  this  source  of  error. 
These  modifications  have  now  been  introduced  into  the  B-2  and  B-3  executive 
modules,  paving  the  way  for  the  subsequent  compilation  of  the  CALIB  program  in 
these  systems. 
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SECTION  4 


SPIRAL  DECAY  EXPERIMENTATION 


During  the  application  at  ADC  of  procedures  developed  under 
this  study  for  spiral  decay  corridor  determination,  six  visual  confirma¬ 
tions  of  decay  were  successfully  scheduled.  Four  of  these  represented 
Russian  satellites  which,  by  virtue  of  their  lower  inclinations,  were 
more  likely  to  decay  over  the  populated  middle  and  equatorial  latitudes. 

Two  of  these  scheduled  visuals  (on  objects  632  and  707)  were  made  during 
the  demonstration  at  ADC  of  the  spiral  decay  capability  in  March,  1964, 
while  a  third  provided  the  first  visuals  on  two  consecutive  revolutions. 

These  three  satellites  are  discussed  in  detail  in  this  section. 

The  visuals  described  above  involved  observers  in  North  America 
and  represent  the  best  documented  observations.  Additional  visuals  were 
scheduled  through  the  1127th  Field  Activities  Group,  Fort  Belvoir,  includ¬ 
ing  object  834,  1  August,  1964,  at  Carnarvon,  Australia;  and  object  796, 

27  May,  1964,  at  Maracaibo,  Venezuela.  The  visuals  on  911  were  made  by 
both  ground  observers  and  SAC  aircraft  commanders,  above  the  Arctic  circle, 
and  were  scheduled  through  NORAD. 

4.1  OBJECTS  632  and  707  DECAY 

To  evaluate  the  performance  of  the  Spiral  Decay  program, 

Milestone  II  version,  in  the  ADC  data  and  computation  environment,  the 
tracking  records  of  two  Russian  Cosmos  Satellites,  Spacetrack  Objects  707 
and  632,  were  reduced  on  a  real-time  basis  in  ADC  centers  at  L.  G.  Hanscom 
Field  and  Colorado  Springs.  Both  exercises  were  successful  beyond  expec¬ 
tations,  resulting  in  visual  confirmation  of  decay  in  each  case.  Equally 
gratifying  was  the  performance  of  the  var iation-of-parameters  method,  coupled 
with  the  analytic  differential  correction  technique  which  (in  the  case  of 
Object  632)  converged  readily  on  tracking  data  extending  into  the  incendiary 
region  of  decay. 

Object  707  was  a  cylindrical  body  of  low  area-to-mass  ratio.  The 
data  reduction  for  this  exercise  was  performed  at  L.  G.  Hanscom  Field  by 
SCAF  and  Aeronutronic  personnel.  The  success  of  this  venture  may  be  measured 
by  results  presented  in  Table  III.  The  low  area/mass  ratio  of  this  object 
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TABLE  III 


CORRELATION  OF  VISUAL  OBSERVATIONS  AND  EPHEMERIS 


SPACETRACK  OBJECT  707  DECAY 
(Based  upon  data  over  3/26  0208Z  -  3/26  2225Z) 


TIME 

(Day:  Time) 

POSITION 

OBSERVER 

COMMENTS 

27: 0128Z 

44°. 2N 

115°. W 

94km 

27: 0129Z 

45°. 7N 

110°.  Ov/ 

89km 

-  - 

27: 0130Z 

47°. ON 

105°. 0W 

82km 

Finley  AF  WX  Station 

Finley ,  N.D. 

"Dim  object  report  moving 

SW  to  NEn 

27: 0131Z 

47°. 9N 

100  .W 

70.0km 

Baudette,  Minn. 

Duluth,  Minn. 

(SAGE  Direction 

Center,  4  observers) 

MRed- orange  object  directly 
overhead;  pieces  falling  off.” 

"Three  objects;  visible  trail  lasting 

5  sec.  behind  first  object;  one  object 
appeared  to  descent  in  local  area." 

27:0132. 54Z 

★ 

Madison,  Wisconsin 
(Moonwatch) 

"Sparking  object,  4°elev  ,  355° 
azimuth" 

*  The  ephemeris  for  object  was  based  upon  a  single  drag  coefficient  implying  structual  integrity. 
The  low  W/CpA  and  confirmed  breakup  would  permit  pieces  to  go  beyond  the  impact  area  corresponding 
to  the  intact  object. 


and  its  shape  suggest  a  collection  of  dense  components  (rocket  motor, 
payload)  joined  by  light  structual  components  and  tankage.  One  would 
expect  the  aerothermodynamic  environment  to  rapidly  strip  away  these 
lighter  structual  components,  creating  a  variety  of  objects  of  differing 
area/mass  ratio,  each  following  an  independent  course.  This  speculation 
is  supported  by  visual  evidence,  which  verified  the  onset  of  high  heating 
rates  (Fargo,  N.D.)>  the  breakup  of  the  satellite  (Baudette ,  Minn.),  and 
the  independent  decay  of  at  least  three  objects  (Duluth,  Minn.).  The 
correlation  between  these  phenomena  and  the  ephemeris  provided  by  the 
Spiral  Decay  program  is  excellent,  as  demonstrated  in  the  Table. 

The  decay  exercise  for  Object  632  was  conducted  at  Colorado 
Springs.  The  area/mass  ratio  of  632  was  also  quite  low,  as  noted  in 
Table  IV.  Nine  fits  were  obtained  over  the  last  72  hours  of  lifetime, 
each  requiring  approximately  10  minutes  of  computer  (Philco  212)  time. 
Highlights  of  each  iteration  are  presented  in  Table iv;  these  include 
(1)  the  fact  that  all  fits  over  the  last  72  hours  (and  probably  greater 
intervals)  defined  a  decay  window  of  approximately  one  hour  and  (2)  no 
iterations  on  any  fit  were  divergent.*  This  latter  property  of  the 
Spiral  Decay  program  is  readily  appreciated  when  it  is  noted  that,  on 
fit  #1  beginning  with  SPADAT  mean  elements,  range-rate  residuals  of 
3  km/sec  and  larger  were  common  on  the  initial  pass  through  the  data. 
Visual  confirmation  of  707  decay  in  its  early  stages  was  provided  by 
a  weather  observer  in  Duluth,  Minn.;  poor  weather  over  the  Eastern 
states  and  mid-Atlantic  precluded  sightings  extending  on  into  the  most 
probable  decay  region  (extending  from  the  Azores  into  South  Africa). 

It  is  highly  probable  that  632  decayed  in  several  fragments 
as  in  the  case  of  707.  In  order  to  evaluate  the  possible  decay  corri¬ 
dor  for  objects  of  various  area/mass  ratios,  a  sequence  of  decay 
trajectories  was  integrated,  beginning  with  the  visual  evidence  of 
decay  over  Duluth  at  0726Z.  In  particular,  a  Moonwatch  observation 
was  reported  about  0900Z  from  the  Northeastern  U.S.  Based  upon  a  fit 
including  tracking  data  during  this  decay  phenomenon,  an  object  with 
an  area/mass  ratio  even  as  low  as  0.0002  m^/kgm  would  not  complete 
another  revolution.  Since  it  is  difficult  to  visualize  a  satellite 
component  with  these  properties,  and  since  both  Moores town  and  Prince 
Albert  radars  gave  negative  reports  on  this  post-decay  revolution, 
confirmation  of  this  report  is  lacking. 


*Since  the  observations  accepted  from  iteration  to  iteration  varied, 
divergence  is  arbitrarily  defined  as  an  increase  in  the  weighted  rms 
of  five  percent  or  more  on  successive  iterations. 
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TABLE  iv 


SPACETRACK  OBJECT  632  DECAY 


PRE -DECAY 

ANALYSIS 

-EUL 

Data  Interval 
(Dav:  Zulu  Time) 

Decay  Time 
(Dav:  Zulu  Time) 

Drag  Parameter 
Meters2  kgm  * 

RMS  Fit 
(Weighted*) 

Divergent 

Iterations 

1 

26:0446  - 

27: 1511 

30:0817 

0.01631 

1.03 

0 

2 

27:1758  - 

28: 1600 

30:0911 

0.01630 

0.70 

0 

3 

27:1753  - 

28:  1732 

30:0856 

0.01647 

0.85 

0 

4 

28:0400  - 

29: 1242 

30:0818 

0.01666 

0.86 

0 

5 

29:0351  - 

29:  1941 

30:0803 

**0.01694 

1.06 

0 

6 

29:0351  - 

29:2113 

*** 

**0.01684 

1.22 

0 
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29:0351  - 

30:0016 

30:0806 

**0.01688 

1  .36 

0 

8 

29:0351  - 

30:0423 

30:0757 

**0.01698 

2.02 

0 

POST-DECAY 

ANALYSIS 

9 

29:1701  - 

30:0727 

30:0747 

**0.01746 
(6  element  fit) 

2.84 

0 

* 

Generally  based 

upon  TIP  bias /weight  recommendations. 

as  revised  by  results  of 

the  object  759 

(1964  March  3-5) 

calibration  effort 

**  Beginning  with  Fit  5,  transitional  drag  effects  were  introduced 

***  Although  solution  was  convergent,  convergence  test  failed  to  meet  stringent  (one  percent)  test  on  successive 
values  of  weighted  rms  before  the  pre-designated  maximum  number  of  iterations  was  exceeded. 


4.2  OBJECT  803  DECAY 


One  of  the  better  documented  decay  exercises  involved  the 
COSMOS  31  (Object  803)  payload,  launched  on  6  June,  1964.  A  TIP 
exercise  conducted  at  SPACETRACK  during  the  week  prior  to  decay 
successfully  scheduled  visual  confirmations  of  decay  during  revolu¬ 
tion  2161,  along  a  corridor  extending  across  the  Northwestern  U.S., 
at  20  October,  0516Z.  Reports  of  unusual  meteoritic  phenomena  on 
the  previous  revolution  2160  prompted  an  appeal  in  Canadian  news 
media  for  detailed  reports. Of  the  resulting  fourteen  reports  sub¬ 
mitted,  a  majority  correlated  with  the  timing  and  direction  of  803, 
leaving  little  doubt  that  some  component  of  803  had  decayed  over 
Ottawa,  Ontario,  at  approximately  2353  local  time  (20  October,  0353Z) . 

Further  analysis  was  conducted  to  determine  answers  to 
the  following  questions: 

(1)  What  mechanism  separated  Object  803  at  least  one 
revolution  before  decay,  and 

(2)  What  was  the  nature  of  the  object  that  decayed 
prematurely? 

The  analysis  presented  here  supports  the  conclusion  that  an  array 
of  solar  cells  was  separated  by  aerodynamic  forces  near  the  ascending 
node,  revolution  2160,  and  decayed  twenty  minutes  later  over  Ottawa. 

Examination  of  the  ephemeris  during  the  last  few  revolutions 
revealed  that  the  perigee  of  the  orbit  was  near  the  ascending  node.  A 
one-minute  interval  ephemeris  was  produced  for  the  last  revolutions  to 
examine  the  altitude  variation  around  the  orbit;  the  following  values 
were  determined: 


Rev 

Altitude  in  Vicinity  of 

Ascending 

Node 

Northern 

Antinode 

Descending 

Node 

Southern 

Antinode 

2159 

141  km 

149  km 

144  km 

157  km 

2160 

133  km 

137  km 

132  km 

142  km 

2161 

119  km 

Decay 

Two  mechanisms  may  be  effective  in  removing  a  component  of 
Object  803  near  the  ascending  node  of  revolution  2160,  aerodynamic  and  ^ 
aerothermodynamic .  The  dynamic  pressure  is  still  quite  low  (~0.025  kg/m  ) 
but  capable  of  inducing  substantial  vibration  loading  on  structures  of 
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poor  aerodynamic  design.  The  aerothermodynamic  effects  at  130  km  border 
on  heating  rates  sufficient  to  induce  structural  failure.  The  combination 
of  these  effects  at  130  km  should  be  sufficient  to  separate  large  surfaces 
with  limited  structure  designed  for  deployment  after  orbit  injection; 
antenna  and  solar  cell  arrays  fall  into  this  category. 

To  evaluate  the  nature  of  the  object  which  decayed  over  Ottawa, 
a  sequence  of  decay  trajectories  was  run  for  a  variety  of  drag  coefficients, 
based  upon  position  and  velocity  derived  at  the  ascending  node,  revolution  2160. 
The  area-to-mass  ratio  leading  to  visual  decay  phenomena  over  Ottawa  (at  110  km) 
was  0.18  m^/kgm,  compared  with  0.01  for  typical  payloads  and  0.03  for  typical 
rocket  bodies.  (A  family  of  these  trajectories  is  plotted  in  Figure  3.) 

The  final  evidence  suggesting  the  actual  identity  of  the  object 
lies  in  the  reported  visual  phenomena.  Some  of  the  reports  follow: 

"...oblong  in  shape,  with  a  long  trail  of  sparks..." 

"...lights  were  a  vivid  yellow  and  seemed  to  be  like  small 
illuminated  windows  in  a  long  row.  Smaller  lights,  like 
sparks,  were  seen  beside  the  main  body .. .approximately 
20  lights  strung  out  in  a  row  resembling  lights  of  a 
passenger  train." 

"...left  a  long  trail  of  illumination  behind  it,  resembling 
the  trail  left  by  an  ascending  rocket  in  a  fireworks 
display . " 

"...a  row  of  lights  down  the  centre ...  looked  like  a  distant 
train  where  the  train1 s  outline  cannot  be  seen." 

"...like  a  stream  of  lights  in  irregular  pattern,  those  lights 
in  front  being  quite  bright .. .Each  light  seemed  to  hold 
its  (relative)  position." 

"...a  continuous  bursting  of  multicoloured  sparks,  which 

appeared  to  be  coming  from  a  common  object  as  they  were 
in  a  straight  line." 

"...lights  were  white  and  surrounded  by  a  yellowish  luminosity 
with  a  long  liminous  tail..." 

Materials  commonly  used  for  optical  enhancement  of  reentry 
objects  include  the  easily-ionized  metallic  oxides  of  such  elements  as 
cesium,  selenium,  indium,  calcium  and  iron.  As  a  corollary,  some  of 
these  materials  are  also  efficient  photoelectric  media  and  are  used  in 
solar  energy  converters.  Thus  the  reported  visual  phenomena  could  be 
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FIGURE  3.  COSMOS  31  (OBJECT  803)  DECAY 


.0109 


explained  by  a  solar-cell  panel,  dispensing  individual  cells  along  the 
higher  trajectory  as  the  bonding  material  failed,  which  acted  as  efficient 
converters  of  heat  to  visual  energy  as  the  lower  altitude  range  was  reached. 
The  area-to-mass  ratio  derived  above  tends  to  support  this  conclusion. 
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APPENDIX  I 


EPHEMERIS  AND  DIFFERENTIAL  CORRECTION 
PROCEDURE 


1.1  VARIATI ON-OF - PARAMETERS  EPHEMERIS  COMPUTATION 

The  formulation  of  the  variation-of-parameters  method  for  ephemeris 
computation  is  presented  in  this  section.*  In  principle,  the  variation-of-para- 
meters  method  formulates  the  equations  of  motion  in  terms  of  dependent  variables 
which,  in  the  absence  of  perturbative  (non- two-body)  forces,  are  constant.  Such 
a  formulation,  used  in  a  differential  correction  procedure,  enhances  the  compu¬ 
tational  efficiency  through  (1)  more  efficient  ephemeris  integration,  particularly 
where  perturbations  are  small,  and  (2)  effective  application  of  differential 
correction  based  upon  analytical  expressions  presented  in  Section  1.2. 

The  variation-of-parameters  method  involves  the  steps  shown  in  Figure 
1-1.  Each  of  these  steps  is  formulated  below;  in  practice,  they  are  programmed 
as  subroutines  which  are  called  by  a  control  program  based  upon  the  numerical 
integration  procedure  used.  Each  pass  through  the  procedure  advances  time  by 
a  designated  amount  At;  repeated  application  will  define  future  position  and 
velocity  from  epoch  elements. 

a .  Derivative  Evaluation 

(1)  Oblateness  (bulge) 

A  fairly  general  Earth  model  has  been  included  in  the  SPIRDEC  and 
CALIB  programs,  inclusive  of  zonal  harmonics  of  the  geopotential  through  the  fifth 
order  and  tesseral  (longitude-dependent)  harmonics  through  the  fourth  order.  The 
most  efficient  procedure  for  computing  these  terms  is  to  define  perturbative 
accelerations  in  the  S  (south),  E  (east)  and  Z  (zenith)  coordinate  system,  as 
follows: 

*B  =  S  gg  +  E  ^  +  Z  gr 

where  gq,  and  g  are  defined,  in  terrhs  of  the  potential  function  $,  as  shown 
in  Table  I-T.  r 


*  See  also  Appendix  VI,  Variation-of-Parameters  for  Low  Eccentricity  Orbits. 


27 


I 


a,  h,  L  at  t 
—  —  o 


FIGURE  1-1.  VARIATION-OF -PARAMETERS  METHOD 
FOR  EPHEMERIS  COMPUTATION 
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Table  1-2  conveys  the  series  to  evaluate  the  associated  Legendre 
functions  and  the  formula  used  to  compute  their  derivatives.  The  derivatives 

always  occur  with  the  factor  Vl  -  U l  in  Table  1-1. 

The  series  expressions  are  not  presently  used  in  the  appropriate  sub¬ 
routine  (MARTINI) .  That  subroutine  conforms  to  the  specifications  in  the 
Milestone  II  Report*,  Appendix  A.  A  new  subroutine  (MONICS)  with  an  extended 
potential  (n^  =  9  and  ^  =  6,  presently)  has  been  validated  but  not  yet 
integrated  into  the  operational  program.  MONICS  uses  the  series,  avoiding  the 
singularities  at  the  equator  and  poles  by  numerical  checks. 

The  values  of  the  sines  and  cosines  of  multiples  of  the  East  longitude 
may  be  found  conveniently  by  the  recursive  relationships: 

sin  n  X£  =  2  sin  (n-1)  XE  cos  X£  -  sin  (n-2)  XE 

cos  n  X£  =  2  cos  (n-1)  X£  cos  X£  -  cos  (n-2)  X£ 

(2)  Drag 

The  perturbations  due  to  atmospheric  drag  are  determined  by 
entering  the  LJDRAG  subroutine , which  determines  the  atmospheric  density 
corresponding  to  the  satellite1 s  altitude  from  a  dynamic  model  atmosphere  and 
translates  this  to  a  corresponding  acceleration.  The  density  model  currently 
employed  is  due  to  Jacchia,  with  approximating  curves  in  transitional  regions 
provided  by  Lockheed.  The  formulation  for  this  density  model  is  given  in 
Appendix  II.  (An  additional  entry  into  this  subroutine  is  used  to  determine 
perigee  density  and  scale  height,  in  order  to  initialize  the  drag  differential 
correction  procedure.) 

To  compute  the  drag  perturbation,  the  height  H  above  the  oblate 
spheroidal  earth  is  first  calculated  (See  Appendix  V,  Section  3): 

H  =  (r-1)  -  |  f2  (u  )4  +  (f  +  |  f2)  U  2 

The  atmospheric  density  corresponding  to  this  altitude  is  then 

determined . 

The  relative  air  velocity  vector,  V,  is  computed  by  assuming  that 
the  atmosphere  is  rotating  at  the  same  angular  velocity  as  the  Earth: 

.  •  •  —  t 

V  x  =  x  4y0  ,  where  8  =  0.058,834,47  radians/ke  min 

U  =  y  -x6 

y 


*  Walters,  L.  G.,  Hilton,  C.  G.,  and  Crossin,  P.  A.,  "Spiral  Decay  and 
Sensor  Calibration  Differential  Correction  Programs",  Aeronutronic 
Publication  No.  U-2559,  31  March  1964. 
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TABLE  1-2.  ASSOCIATED  LEGENDRE  FUNCTIONS  AND  THEIR  DERIVATIVES 


Legendre  Polynomials:  m  =  0 

P  (U  )  =  P  n  (U  ) 
n  z  n,0  z 


Associated  Legendre  Functions: 
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Derivatives: 
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The  drag  acceleration  then  is 


where  B 


/m  is 


Pv  B  u 

the  ballistic 


parameter . 


(3)  Radiation  Pressure 

The  magnitude  of  the  force  acting  on  a  satellite  due  to  direct 
solar  radiation  pressure  is 


F 

© 


P  A 
© 


where  A  is  the  effective  cross-sectional  area  of  the  satellite  to 

radiation  pressure 


P  is  the  solar  radiation  pressure  in  the  vicinity  of  the 

Earth,  assumed  constant  at  4.5  x  10  5 

2 

cm 


7  is  a  factor  depending  on  the  reflecting  characteristics  of 

the  satellite's  surface;  if  the  incident  energy  is  reflected 
specularly  or  is  absorbed  and  re-emitted  isotropically, 

7  =  1 


By  neglecting  the  solar  parallax,  which  is  only  11  seconds  of  arc  at  1000  miles 
altitude, 


-F  L 


0  ~ © 


where  L 


is  a  unit  vector  directed  to  the  Sun. 


The  points  at  which  a  satellite  enters  and  leaves  the  Earth's 
shadow  are  continuously  changing  due  to  the  apparent  motion  of  the  Sun 
and  the  perturbations  on  the  orbit.  The  shadow  effect  can  be  handled 
very  easily,  however,  in  the  special  perturbations  case.  A  check  is  made 
at  each  integration  step  to  determine  whether  the  satellite  is  illuminated 
or  in  the  Earth's  shadow.  If  the  satellite  is  in  the  shadow,  the  radiation 
force  does  not  act. 

Consider  the  plane  which  passes  through  the  centers  of  the  Earth, 
Sun,  and  satellite  as  shown  in  Figure  1-2.  The  Earth's  shadow  is  assumed 
to  be  cylindrical.  The  dot  product  of  I^and  r  gives 

.  L*  r 

cos  w  =  — 0  — 

r 

Clearly,  if  cos  if  is  positive,  the  satellite  will  be  illuminated  by  the 
Sun.  If  cos  i/  is  negative,  use  is  made  of  the  triangle  CES  and  the 
angle  V,  obtained  from  triangle  TES: 

.  *n  1 

sin  V  =  — 
r 

Consideration  of  Figure  1-2  reveals  that  when  (ilr  +  T)  180°,  the  satellite 
is  still  illuminated.  However,  when  (i|r  +  T)  >  180°,  the  satellite  will  be 
in  the  Earth's  shadow.  In  this  case,  the  radiation  force  does  not  act. 

The  test  can  readily  be  made  on  the  sign  of  the  quantity 

sin  (if  +  TO  =  sin  if  cos  V  +  sin  V  cos  if 
where  sin  if  and  cos  Tl  are  always  positive. 
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Figure  1-2  Position  of  Satellite  With  Respect  to  Earth’s  Shadow 
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The  ecceleretlon  due  to  direct  eoler  redletlon  pressure. 

Is  computed  as  follows.  Initially,  the  teeen  longitude 

of  the  sun  et  the  beginning  of  the  current  peer  is  found  by  the  systen 
subroutine  TLC 
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If  sin  VP  +?))  >  0,  the  satellite  is  Illuminated. 
If  the  satellite  Is  Illuminated,  calculate 
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b .  Calculation  of  Variation  of  Parameters 


The  total  perturbative  acceleration  can  now  be  used  to  determine  the 
perturbative  derivatives  of  the  parameters,  as  follows: 
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The  derivatives,  as  used  in  the  Runge-Kutta  routine,  are: 
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c .  Numerical  Integration  Procedure 

The  expression  given  above  represent  seven  first-order  differential 
equations,  with  time  as  the  dependent  variable.  In  the  Milestone  I  and  II 
versions  of  the  SPIRDEC  and  CALIB  programs,  these  are  integrated  by  the  Runge- 
Kutta  procedure,  utilizing  an  optional  Simpson's  rule  test  for  truncation 
error  to  define  an  optimum  variable  time  interval. 


d .  Calculation  of  Position  and  Velocity 

At  each  time  point  where  derivatives  are  to  be  evaluated,  or  where 
observations  are  to  be  represented,  the  space  position  and  velocity  are 
determined  from  the  parameters  a,  h  and  L  by  the  following  procedure: 
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Define  the  nodal  vectors  M  and  N 
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Kepler’s  equation  is  solved  by  iteration  using  the  Newton-Raphson 
thod  with  an  initial  guess  for  (E+  u)  )  of  U  [  i.e.,  (E+  ^  »  U  ] 
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where 


e  cos  E.  *  a  cos  (E+u)).  +  a  sin  (E+u)  ) 
i  xN  1  yN  1 

e  sin  E.  «=  a  sin  (E+u)).  -  a  XT  cos  (E+  a))  . 

1  xN  l  yN  i 

The  iteration  is  concluded  when 

(E+cu)i+1  -  (E+u)). 

(If,  after  50  iterations,  the  criterion  is  not  met,  the  run  is  terminated.  A 
comment  to  this  effect  is  written  on  the  output  tape.) 

After  Kepler* s  equation  is  solved,  the  calculations  continue: 
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U  *  N  cos  u  +  M  sin  u 

V  *  -  N  sin  u  +  M  cos  u 

r  =  r  U 

t  ■  r  U  +  rv  V 

This  completes  the  calculation  of  position  and  velocity  from  the  M,  N  parameters 
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1.2  WEIGHTED  DIFFERENTIAL  CORRECTION 


The  differential  correction  procedure  relates  observation  residuals 
to  incremental  changes  in  the  orbit  parameters.  Raw  observation  data,  corrected 
for  known  systematic  (bias)  errors,  are  compared  against  their  representations, 
based  upon  the  computed  orbit.  Their  differences  (residuals)  are  related  by 
scalar  differential  expressions  (equations  of  condition)  to  incremental  changes 
in  the  parameters  defining  the  computed  orbit.  This  ordinarily  heavily  over¬ 
determined  system  of  equations  is  solved  for  parameter  changes  in  the  sense  of 
weighted  least  squares,  utilizing  weights  reflecting  the  statistical  confidence 
in  the  corrected  data,  as  determined  from  sensor  calibration  efforts. 

.  a .  Representation  of  Observations 

Given  the  sensor  latitude  $,  east  longitudeXg  and  sea-level  height 
H,  the  following  procedure  computes  representations  of  range  p  ,  range-rate 
pc  ,  and  the  direction  cosines  of  the  unit  vector  L  directed  from  the  sensor 
to  the  satellite  for  time  t,  in  minutes  since  epoch: 

Compute  local  sidereal  time  0  at  time  t: 

9  r  +  0  .  +  0.004,375,269,1  t  (Mod  2  n  ) 

Cl  to 

Compute  sensor  location  vector,  R: 

X  -  -  £  (1-  e  ^  sin^0)  ^  +  H  1  cos  0  cos  0 

Y  -  -  [  (1-  sin^0)  ^  +  H  J  cos  0  sin  0 

Z  «  -  (1-  c  2)  (1-  e  ^  sin^0)  ^  +  H  J  sin  0 

Compute  the  slant  range,  p^: 

*  £  +  R,  where  £  is  the  satellite’s  computed  position  at  time  t. 
Pc  =  <£c  '  £c)* 

Compute  unit  vector  from  the  station  to  the  satellite  in  the  equatorial  co¬ 
ordinate  system: 
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Compute  range  rate,  p 


Pc"i  +  ^ 


,  from 

where  the  components  of  R  are  given  by 


/ 


X  »  -Y  9,  0  =  0.058,834,47 

R  /  Y  -  X  0 


Z  -  0 

L 

then  p  c  *  L  *  £  Q  ^  \  (x  +  X)  +  Ly  (y  +  Y)  +  L^z , 

b .  Definition  of  Residuals 

If  range  is  observed,  the  residual  is 


»!  -  P  -  Pc 

If  azimuth,  A,  and  elevation,  h,  are  observed,  the  residuals  are, 
respectively , 

r2  *  p  c  ?  •  (L  -  Lc) 

R3  -  p  c  D  •  (L  -  Lc) 


where 


A  =  A,S  +  A,E  +  A,Z 

—  xh  —  yh  —  zh  — 

D  =  D  ,  S  +  D,  E  +  D,  Z 

—  xh  —  yh  —  zh  — 


L  =  L,S  +  L,E  +  L,Z 
—  xh  —  yh  —  zh  — 
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The  $>  E,  Z  unit  vector  system  and  the  horizon  oriented  unit  vector 

system  are  defined  by: 


S  *  sin  i  cos  0 
x 


S  <  S  s  sin  sin  0 

~  '  y 


,  S  =*  -  cos  i 
z 


M 


"xh 


yh 


-  cos  A  cos  h 
sin  A  cos  h 


L  ,  =  sin  h 
zh 


f E* 


-  sin  0 


E  {  E  ■  cos  0 

y 


/  ^ 

A  ,  “  sin  A 
xh 


A  ,  =  cos  A 
yh 


IE  »  0 
z 

Z  ■  COS  i>  cos  0 

X 

Z  \  Z  *  cos  6  sin  0 

y 

*  Z  •  sin 
z 

If  right  ascension,  a 
residuals  are: 


and  declination ,  6 


-  0 

*  cos  A  sin  h 

»  -  sin  A  sin  h 

»  cos  h 

are  observed,  the 


*  P  c  —  ~ 

r5  "  P  c  D  *  (L  -  Lc) 


where 


L 


L  *  cos  6  cos 

X 


L  ■  cos  6  sin 

y 


L  *  sin  6 
z 


a 

a 
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sin  a 


A 


cos  a 

0 


D  =  -  sin  6  cos  a 
x 


D  <  D  ■  -  sin  a  sin  6 

-  \  y 

,  D  =  cos  6 

'  Z 

See  Figure  1-2  for  these  vector  relationships 

If  range  rate,  p  ,  is  observed, 

R,  -pAp-(p-p)  p 
o  c  c  c 

c .  Differential  Correction 


The  residuals  in  the  observations  are  related  to  the  changes  in 
the  elements  through  first-order  scalar  differential  expressions  of  the  form 


\ 

Anc  /  \ 

/  \ 

) 

c 

1  - 2  +  c  . 

A  a  +  1 

c  .  ' 

A  a  + 

An. 
n  i 

1  \  xnJ 

xn 

i  °  i 

A  a 

i  ynl 

o 

c 

+  C 


/ 

AO  +  ;  C 


ni  ° 


A  i  + 
o 


AB 
B  It 


AJI 

B 


where  the  form  of  the  coefficients,  depend  on  the  observation  type,  time 

of  observation,  and  the  observation  weights.  These  coefficients,  which  have 
been  developed  by  means  of  first-order  partials,  are  functions  of  the  orbit 
elements  and  computed  observations. 


The  coefficients  are  computed  by  first  establishing  the  R  and  U 
coefficients  at  time  t. 
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Celestial  Spher 
Centered  At 
Observer 


FIGURE  1-3 .  VECTOR  RELATIONSHIPS 


R  =  (—  )  e  sin  E 
u  r 


R  =  -f  r  +  (U  -  U  )  R 
n  3  o  u 


RxN  -  <f >  '  “,N  cos  (E  +  u)  )  J 
RyN  '  >  t  ay»  -  Sl"  (E  +  "  )  i 


U  =  (-  )  V  1  -  e2 


u  r 


U  =  (U  -  U  )  U 
n  o  u 


U  M  =  -  <(!+—)  sin  (E  +  cu  )  + 


xN  r 


a  e  sin  E 
xN 


’  e2  -  (1  +  •[  1-e2 


(1 


(1  -f  V  1-e  )  e  cos  E 

+  Vl-e2)2  Vl-e2 


yN 


1  +  Vl-e 


U  „  -  -  <-  (1  +  -  )  cos  (E  +  u)  )  + 


a 

"yN  r 


a 

r 


a  XT  e  sin  E 
yN 


e2  -  (i  +  yr 2 


e  )  e  cos  E 


xN 


(1  +  V  l-e2)2  Vl-e2 

(1)  Coefficients  for  slant  range  observation 


1  +  Vl-e 


When  the  slant  range  is  observed,  the  residual  and  coefficients 
are  given  by  the  expressions: 


R  ■  R.  ct 
1  P 


-1 
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c 


.  -  r  L  -UR+L  •  V  U  1  a 

A  n  — c  —  n  — c  —  n  J  | 


-1 


n 


A  a. 


xN 


yN 


A  U 


A  0 


A  i 


[  L  '  URh  +  L  •  V  U  H  1  0~ 1 

[  — c  —  xN  — c  —  xN  J  p 

[L  -  URm  +  L  •VUMlcr"1 

l  -c  —  yN  — c  -  yN  J  p 

[l  -ur+l  -  vuIct'1 

[*. 


V  r  cos  i  -  Lc  ’  W  r  sin  i  cos 


“] 


-1 


W  r  sin 


“] 


The  expressions  for  ^  are  presented  in  Appendix  iv. 

(2)  Coefficients  for  azimuth  and  elevation  observations 

When  azimuth  A  is  observed,  then  R  »  R  (  p  q  » )  and  the 

^  4  C  A 

coefficients  are  obtained  by  replacing  L  by  A  and  by  p  q  a  where  a  ^ 
is  the  standard  deviation  in  the  azimuth  measurement. 

For  elevation  observations,  R  -  R.  (  p  a  ,  )  *  and  the  corre- 

3  c  n 

sponding  coefficients  are  obtained  by  replacing  L  with  B\  is  replaced  by 

p  c  cj  h>  where  a  K  is  the  standard  deviation  in  elevation  angle. 


puted  from: 


(3)  Coefficients  for  range-rate  observations 

For  range  rate  observations,  preliminary  coefficients  are  com- 


u 


<-£)*  (?) 3  ( 


e  cos  E  -  e  ) 
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*N 

-  |  +  (u 

R  „ 
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xN  a 


<—  )*  (7) 3  V  1  -  e2 


U 


yN 


ajs 

n 


<—  (r)3  7  1  -  e2 


2\ 
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1  + 

—  i 

flP/  J 

- 

1 

f 

2  ^  1 
y» 

sin  (E  +  a)  )  -  a^N  1 

1  + 

\ 

L 

ap 

•* 

are  then  computed  as 

)  -  p  R  1  +  p  *  U 
iy  KcnJ  —  c  — 

R 

n 

taxN  1  "c 


+  V  V  [pc(Un  V  S  Un  [  °p' 

Uo-at  »c  <Rx»-'' V  -  PcR*N]+  £c‘  -  R> 


xN 


+  £c  •  2C  Pc  (“*H  +  V  U*N  •  Pc  “xN  5  +  £c‘  -  U 


xN 
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The  coefficient  C  is  given  in  Appendix  IV. 

A  B 
~5~ 

(4)  Differential  Correction  Solution 


Let 


£ l  °ij  A  j  "  Rij 


represent  all  such  equations  of  condition,  where  C  are  the  coefficients, 

R  are  the  accepted  observation  residuals,  A.  areJthe  corrections  to  the 
orbital  elements,  and  N  is  the  number  of  elements  to  be  corrected  and  is  the 
number  of  accepted  observation  residuals.  The  resulting  matrix  equation  is 
solved  to  give  the  correction,  A  .in  a  least  square  sense,  to  the  orbital 
elements  at  time,  t  .  These  corrections  are  applied  as  follows  (primes 
denote  corrected  elements): 
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xN 
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yNo  yNo  yN 
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o  o 
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! 
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°  °  Z 
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o  o  o  z 

\ 

Following  the  above  calculation  of  the  corrected  elements,  another 
representation  of  the  observations  is  performed,  on  the  basis  of  the  corrected 
elements,  and  another  set  of  residuals  is  formed  by  using  the  same  input 
observations.  The  weighted  RMS  values  of  sets  of  consecutive  residuals  are 
compared  to  define  convergence  of  the  computational  process. 
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APPENDIX  II 


JACCHIA-NICOLET 
DYNAMIC  ATMOSPHERE  MODEL 


The  physical  properties  of  the  outer  atmosphere  are  governed  by 
two  relations  between  pressure  P,and  density,  p.  The  first  is  the  equation 
of  hydrostatic  equilibrium 

dP 

dh  "gp 

where  g  is  the  acceleration  of  gravity  at  height  h.  The  second  is  the  ideal 
gas  law. 


MP  =  kpT 

Here  k  is  the  Boltzmann  constant,  but  M,  the  mean  molecular  mass,  and  T,  the 
temperature,  are  additional  variables.  Fortunately,  the  temperature  is 
practically  independent  of  h  above  a  level  called  the  thermopause.  The 
variations  of  mean  molecular  mass,  however,  must  be  assumed  according  to  some 
theory.  This  theory  concerns  itself  with  the  dissociation  of  molecules  and 
the  diffusion  of  atmospheric  constituents.  The  resulting  atmospheric  model, 
obtained  by  integrating  the  above  equations  is  a  static  atmosphere.  It 
depends  also  on  the  boundary  conditions  assumed  at  the  lower  boundary.  Upon 
this  static  model,  variations  due  to  solar  activity  can  be  effected. 
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By  employing  Nicolet's  model  for  density  variations  with  height 
and  temperature,  Jacchia^  has  demonstrated  excellent  correlation  between  solar, 
geomagnetic  and  geographic  phenomena  and  exospheric  temperature,  the  latter 
derived  from  satellite  accelerations  using  Nicolet’s  static  model.  By  reversing 
this  procedure,  temperature  may  be  derived  using  Jacchia's  formulas,  and  the 
density  from  Nicolet's  model;  this  procedure  has  been  implemented  in  the  Spiral 
Decay  Program  to  provide  a  flexible  dynamic  atmosphere  model. 

II ,1  RELATION  OF  SOLAR,  GEOMAGNETIC  AND  GEOGRAPHIC  PHENOMENA  TO  TEMPERATURE. 

Temperature  is  computed  from  Jacchiafs  expressions,  beginning 
with  the  input  quantities. 

F1(J  10.7  cm  flux 

Fio  F^  averaged  over  three  solar  rotations 

Ap  Geomagnetic  index 

The  following  are  provided  by  the  program: 

D  Day  number 

Oq,  ^0  Position  of  the  sun  (right  ascension  and  declination) 

0  Sidereal  time  at  vehicle 

0  Geocentric  latitude  of  vehicle 

h  Height  above  sea  level 


Nicolet,  M.,  "Density  of  the  Heterosphere  related  to  Temperature"  SAO  Special 
Report  75,  19  September  1961.  (SAO  has  furnished  an  improved  tabulation  of 
this  model.) 

o 

Jacchia,  L.  G.,  "The  Temperature  above  the  Thermopause"  SAO  Special  Report  150, 
22  April  1964. 
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Temperature  is  computed  as  follows: 

Compute  the  average  night  time  minimum  temperature  (°K) 

T  =  974°  +  4 ,°2  -  150)  +  0.°004  (Fin  -  150)2 

o  10  10 

Compute  the  night  time  minimum  for  the  day 

To  -  \  +  1-°9  (F10  -  ?10> 

Add  the  semiannual  term 

T  =  t'  +  [  0.°39  +  0.°15  sin  2"  f  D~~) 
o  o  [  V  365  ;J 

x  sin  4^ 

Correct  for  latitude  and  elongation  from  sun  (see  figure  II-l) 
Define:  2  T)  =  (0 '  -  60) 

2  S  =  (0'  +  a0) 

the  maximum  daytime  temperature  at  latitude  0/  is 

T  =  T  (1  +  R  cos'"  71) 

D  O 

and  the  minimum  nighttime  temperature  at  latitude  07is 
TN  =  To  (1  +  R  sin"?) 

where  Tq  is  the  global  minimum  in  the  exosphere. 

Empirical  values  for  R  and  m,  derived  by  Jacchia  and  Slowey,  are  0.3 
and  2.5  respectively. 


I D-60\ 

I  365  I 


Hg  =  d  -  otg  (radians) 

T  £  H0  -  7  +  0.21  sin  (Hn  +  y)  -  tt  <  t  <  +  tt 
0  4  0  4 

T'  -  To  -  0.3  Tq<J  sin  2*5|  j  1  -  cos  2,5  £  j+  cos  2,5T|  cos 

Finally  the  geomagnetic  effect  is  added 
T  =  T'  +  l.°2  A 

P 


2.5  T 
2 
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Cone  Containing  Vehicle  Radius  Vector 


FIGURE  II- 1.  RELATION  OF  VEHICLE  LATITUDE  TO  SUBSOLAR  BULGE 


II. 2  APPLICATION  OF  NICOLET'S  STEADY-STATE  MODEL 


Given  the  altitude  and  temperature,  density  is  determined  by  inter¬ 
polation  in  Nicolet's  steady-state  model.  This  model,  tabulated  for  10  km 
altitude  and  50°  temperature  increments,  involves  almost  3000  storage  locations. 
In  order  to  reduce  this  number  to  a  reasonable  value,  series  representations 
were  evaluated.  In  particular,  for  each  altitude,  a  third  order  polynomial 
provides  a  satisfactory  representation  for  the  variations  of  log  density  with 
temperature,  i.e. 


r 

/t  rri\  /t  1  1  AA\  *  T"  1100 

+  B 
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-l-  O 

T-1100 

p  (h ,T) -  log  p  (h , 1100)  -  A  ^00 

400 

T-  L 

400 

w  - 

4- 

The  four  values  needed  to  determine  log  density  from  temperature  --  log  p  (h,1100), 

A,  B,  C  --  for  each  altitude  increment  are  compacted,  two  words  to 
each  cell,  into  the  48  bit  word  structure  of  the  Philco  2000,  requiring  less 
than  300  core  locations.  Below  120  km,  the  limit  of  the  Nicolet  II  tables, 
values  of  log  p  derived  from  the  COESA  1962  model  are  used;  the  densities  at 
100  and  110  km  have  been  adjusted  slightly  to  assure  a  smooth  transition.  At  120  km, 
the  SAO  density  has  also  been  decreased  slightly.  This  is  in  the  direction 
indicated  by  Slowey's  reduction  of  Explorer  17  data.3 

Given  the  temperature,  four  values  of  density  are  derived  for  two 
altitude  increments  above  and  below  the  actual  satellite  altitude.  Density 
and  scale-height  for  the  satellite  altitude  are  derived  from  these  points. 

This  model  follows  Nicolet's  selection  of  molecular  temperature  as 
the  essential  parameter.  Any  refinements  in  the  relationship  between  this 
temperature  and  geophysical  parameters  can  easily  be  incorporated  into  the 
model  in  the  future.  For  instance,  the  indication^  that  the  geomagnetic  index, 

K  ,  is  better  correlated  with  atmospheric  fluctuations  than  is  A  may  result  in 
a^different  "geomagnetic  effect"  formula.  ^ 

II. 3  INTERPOLATION  FORMULAS 

Given  the  temperature,  four  consecutive  values  of  log  p  are  derived, 
two  for  altitudes  higher  than  the  satellite  altitude  h,  and  two  below.  (Near 
the  ends  of  the  table,  more  values  will  be  used  on  one  side  than  the  other). 

These  altitudes  are  designated  by  the  subscripts  -1,  0>  1>  and  2.  Then, 

defining  s  =  0.1  (h-hg  )  for  10  km  tabular  intervals  in  altitude, 


^Slowey,  J.  "Atmospheric  Densities  and  Temperatures  from  the  Drag  Analysis  of 
the  Explorer  17  Satellite"  SAO  Special  Report  157,  1  July  1964. 
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log  0  (h,T)  »  log  p  (h_1>  T) 

+  (s+l)t,-l)<>-2?  log  p  T) 

-  log  0  (hj,  T) 

+  istlHsKs^l  logp(h2>T) 


The  density  scale  height,  H  ,  is  given  in  km.  by 

P 

l  A  log  p 

H  ’  A  u6 

p  An 


l°ge10 

10 


3s  -  6s  4*  2 


log  p  (h_1,  T) 


+  3s  -  4s  -  1 


log  p  (hQ,  T) 


-  3s  -  2s  -  2 


log  p  (h^  T) 


3s2  -  1 

-6  -  log  p  (h2,  T) 


II. 4  DENSITY  ABOVE  1000  km 


Above  1000  km,  the  density  is  extrapolated  from  tabular  values  at 
1000  km  by  assuming  an  isothermal  atmosphere  with  constant  scale  height  as  de¬ 
termined  at  1000  km,  i.e. 

log  p  (h,  T)  =  log  p  (1000,  T)  -  I^OOOTt)  log  6 

Thus,  the  normal  interpolation  procedure  is  used  to  derive  H  (1000,  T)  from  four 
values  of  log  p  (h,  T)  for  h  =  970,  980,  990,  and  1000  km. 
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Since  s  =  2, 

log  10  e 

H(1000,T) 


2 

Z  a4  lo8  P  0*.  »  T) 


where 


-1 


970  km 
980 
990 
1000 


-1 


+  1/30 

-  3/20 
+  3/10 

-  11/60 


Thus,  above  1000  km 

2 

log  p  (h,  T)  =  log  p  (1000,  T)  -  (h  -  1000)  Z 

i=  - 


a±  log  p  (h±,  T) 
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APPENDIX  III 

SCALAR  DIFFERENTIAL  EXPRESSIONS  FOR 
ORBIT  PARAMETERS 


An  orbit  determination  is  a  procedure  for  translating  observations 
into  a  theory  of  the  orbit  suitable  for  producing  a  prediction  or  ephemeris. 
The  orbit  determination  method  employed  herein  is  based  upon  differential 
correction  preceded  by  a  "representation”  or  calculation  of  the  observations 
from  a  preliminary  or  approximate  theory.  The  representation  is  accomplished 
by  special  perturbations  or  numerical  integration  of  selected  equations  that 
define  the  motion  or  departures  from  an  initial  reference  orbit.  In  the 
interests  of  computational  efficiency,  two  departures  from  the  engineering 
practice  of  differential  orbit  correction  have  been  made: 

(1)  The  linear  differential  correction  formulae  are 
obtained  analytically  by  differentiating  relation¬ 
ships  between  the  observations  and  the  orbit,  and 

(2)  The  representation  of  the  observations  from  a  pre¬ 
liminary  orbit,  required  for  the  computation  of 
observation  residuals,  is  carried  out  by  the 
variation-of-parameters  method,  thus  avoiding  the 
time-consuming  aspects  of  the  numerical  integration 
of  the  total  acceleration  into  an  orbit. 
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The  resulting  method  has  demonstrated  outstanding  computational 
efficiency,  as  evidenced  in  the  experimentation  described  in  a  subsequent 
volume,  and  ability  to  cope  with  a  wide  variety  of  observation  patterns. 

This  section  deals  with  the  choice  of  parameters  utilized  in 
the  orbit  theory  and  presents  a  derivation  of  the  linear  correction  form¬ 
ulae  for  translating  topocentric  observation  residuals  into  corrections 
of  the  orbit  parameters  utilized  in  the  representation  of  the  observa¬ 
tions. 


III.  1  SELECTION  OF  ORBIT  PARAMETERS 

The  degree  of  success  of  a  differential  correction  depends 
largely  upon  the  perspicacity  of  the  choice  of  the  parameters  that  rep¬ 
resent  the  orbit  in  the  correction  process.  Any  choice  of  parameters 
which  leads  to  nearly  linear  relationships  between  cause  (errors  in  the 
parameters)  and  effect  (topocentric  observation  residuals)  over  the 
expected  range  of  the  values  of  the  parameters  forms  the  basis  for  an 
adequate  orbit  theory.  Additional  considerations  of  computational 
simplicity  and  compatibility  with  other  aspects  of  the  representation 
or  correction  procedures  serve  further  to  restrict  appropriate  choices 
of  the  parameters.  A  conventional  choice,  including  eccentricity,  e, 
and  argument  of  perigee,  CD,  is  unsatisfactory  because  of  the  inevitable 
tendency  toward  lower  eccentricities  of  geocentric  satellites  interacting 
with  the  upper  atmosphere.  As  the  eccentricity  approaches  zero,  prior  to 
decay,  the  argument  of  perigee  becomes  more  poorly  defined  and  loses  its 
value  as  a  reference  direction  from  which  the  position  (anomaly)  of  the 
object  is  measured.  Alternatively,  a  choice  of  parameters  involving 
nodal  longitude,  ,  and  inclination,  i,  leads  to  a  singularity  as  i 
approaches  zero,  for  Q,  becomes  poorly  defined.  The  latter  consideration 
is  not  important  in  the  present  application  but  must  be  considered  prior 
to  the  appearance  of  equatorial  satellites. 

The  choice  of  parameters  utilized  in  the  development  of  this 
theory  involves  quantities  identified  with  a  coordinate  system  based 
upon  orthogonal  unit  vectors,  N,  M,  W  (Figure  III- l)directed  respectively 
to  the  ascending  node,  the  northern  ant inode ,  and  normal  to  the  orbit  plane. 

The  six  parameters  selected  are  the  following: 

L  mean  longitude  of  the  satellite  measured  in 

two  planes  from  the  vernal  equinox  (i.e. 

n  +  CD  +  M  ) 
o 

a  „  =  e  cos  OQ 
xN 


a  „  =  e  sin  0) 
yN 
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NORTH 

CELESTIAL  POLE 


FIG.  Ill- 1. PROJECTION  OF  ORBIT  ON  CELESTIAL  SPHERE,  WITH 
ORIENTATION  UNIT  VECTORS  AND  ANGLES  DISPLAYED 
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a  semi-major  axis  of  Ihe  osculating  orbit 

longitude  of  the  ascending  node 
i  inclination  of  the  orbit  from  the  equator  plane. 


The  choice  of  the  M,  N  and  W  coordinate  system  and  of  Cl  and  i  as  parameters 
are  subject  to  revision  when  equatorial  orbits  are  considered.  These  param¬ 
eters  are  wholly  adequate,  however,  for  the  satellites  which  are  the  subjects 
of  the  present  study. 

Another  argument  governing  the  choice  of  orbit  parameters  is 
the  compatibility  with  the  ephemeris  integration  program.  In  the  variation 
of  parameters  procedure  employed  in  this  study,  the  derivatives  of  LQ,  £ 
and  the  angular  momentum  h  =  Jp  are  determined  from  the  perturbative 
(non  two-body)  force  field;  their  integrals  then  do  not  involve  the  domi¬ 
nant  central-mass  gravitational  term  and  an  extremely  efficient  ephemeris 
program  results.  The  duplication  of  parameters  between  the  correction 
formulae  and  the  ephemeris  integration  formulae  leads  to  minimum  of  trans¬ 
lation  steps  between  the  two. 


III. 2.  THE  DIFFERENTIAL  RELATIONSHIPS 

This  presentation  is  logically  divided  into  three  parts.  The 
first  is  a  derivation  of  the  differential  expressions  relating  position 
and  orbit  parameters  in  the  orbit  plane  and  involves  the  properties  of 
a  two-body  orbit.  The  second  extends  the  differentiation  to  the  orienta¬ 
tion  of  the  orbit  plane  (specified  by  Cl  and  i  in  this  formulation).  The 
third  step  brings  in  the  observer  and  leads  to  the  final  differential 
correction  formulae  employed  in  this  study. 

a .  Differential  Expressions  for  Position  Components  in  the  Orbit  Plane 

This  section  is  concerned  with  the  derivation  of  the  differential 
relationships  between  the  adopted  parameters  and  the  position  components  in 
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the  orbit  plane.  These  derivations  are  conveniently  separable  from  the 
broader  problem  in  that  these  quantities  depend  only  on  four  of  the  six 
parameters,  i.e.,  a,  UQ  =MQ  +  0)  =  L0  -ft,  ax^  =  e  cosq),  ayN  *  a  sin  q) 
and  the  analysis  involves  only  scalar  two-body  formulae.  The  remaining 
differential  relationships  are  derived  in  section  III, 2b  from  the  expres¬ 
sions  for  the  vectors  defining  the  orientation  of  the  orbit  plane. 

The  whole  of  the  analysis  will  employ  the  N,  M  and  W  unit 
vectors  to  specify  the  orientation  of  the  orbit  plane  and  the  reference 
directions  in  this  plane.  The  N  and  M  vectors  lie  in  the  orbit  plane 
with  N  directed  to  the  node  as  shown  in  Figure  III-l.The  sense  of  W  is 
determined  by  that  of  the  angular  momentum,  northerly  for  direct  (east¬ 
ward)  motion.  Position  in  this  plane  is  denoted  by  the  components  x ^ 
and  y.^  defined  by: 


=  r  cos  u  =  N-r  (1) 

y^  =  r  sin  u  =  M*r  (2) 


where  u,  the  "argument  of  the  latitude,"  is  given  by 


u  =  v  +  a) 


(3) 


These  angles  and  other  quantities  used  in  this  section  are  defined  in 
Figure  IH-1. Other  two-body  expressions  used  in  this  development  include 

n  =  k'Vjiia~  (4) 


U  =  U  +  n  (t  -  t  ) 
o  v  o 

E  —  e  sin  E  =  M  =  U  -  CD 

r  =  a  (1  —  e  cos  E) 
r  cos  v  =  a  (cos  E  —  e) 

r  sin  v  =  a  sin  E 


(5) 

(6) 

(7) 

(8) 
(9) 
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Vlf 

Reference  to  the  glossary  will  clarify  the  meanings  of  these  quantities. 

The  differential  expressions  for  Xjg  and  y^  are  obtained  by 
differentiating  (1)  and  (2). 


Axn  ■  Ar  cos  u  —  rA  u  sin  u  (10) 

Ay^  =  Ar  sin  u  4-  rA  u  cos  u  (11) 

The  determination  of  Ar  follows  from  (7).  The  determination  of  rAu 
is  made  in  terms  of  its  components  rAv  and  rAcD,  the  first  of  these 

from  the  derivatives  of  (8)  and  (9)  and  the  latter  from  the  derivatives 

of  the  parameters  e  cos  cjd  and  e  sin  CD.  These  parameters  will  be  denoted 
axN  anci  ayN>  respectively,  in  subsequent  developments.  They  are  the 
components,  referred  to  axes  determined  by  N  and  M,  of  a  vector  £  directed 
to  perigee  with  magnitude  e. 

The  procedure  outlined  above  is  relatively  straightforward  and 

will  be  summarized  in  this  section.  Starting  with  the  parameters  a^ 

and  a  „  : 
yN 


Aa  =  A(e  cos  (D)  =  Ae  cos  O)  —  e  Ad)  sin  CD 
xN 


Aa  =  A(e  sin  CD)  =  Ae  sin  CD  4  e  Acd  cos  CD 


Expressions  for  Ae  and  eAco  follow  directly: 


Ae  =  Aa  _ _  cos  CD  4  A  a  sin  CD 
xN  yN 


e  Acd  =  —A a  sin  CD  4  Aa  VT  cos  CD 
xN  yN 


Next,  the  differential  of  (7)  yields  Ar: 


r  A  a 

Ar  =  - -  aeAE  sin  E  —  aAe  cos  E 

a 


The  term  AE  follows  from  (4),  (5),  and  (6),  for 


*  Appendix  VII 


(12) 

(13) 


(14) 
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A  n  _  _  3_  A  a 
n  ~  2  a 


AU  -  AUo  +  (t  -  tD)  An  -  AU,  -|(U-Uo)  Ai 


and  finally 

AE  =  -  AU  -  -  AO)  -  ^  (U  “  U  )  -  —  +  -  sin  E  A  e  . 
ror  2  o  rar 


Substitution  for  AE  in  (14)  leads  to 


Ar  =  [  r  -  ^  (U  -  U  )  —  e  sin  E  1  +  AU  \  —  e  sin  E  ] 

a  L  2  err  J  olr  J 


+  Aod  [—  --  e  sin  E  ]  +  Ae  [  —  a  cos  E  +  ^-  e  sin^  E  ] 


Further  substitution  for  Aod  and  Ae  from  (12)  and  (13)  leads  to  the 
form 


Ar  =  R  AU  +  R  —  +  R  „  Aa  „  +  R  „  Aa  M 
u  o  a  a  xN  xN  yN  yN 


(15) 


The  R's  are  the  partial  derivatives  of  r  with  respect  to  the  indicated 
parameters  and  are  tabulated  below: 


R  =  —  e  sin  E 
u  r 


R  =  r-  -(U-U)R 
a  2  v  o'  u 


R  M  =  —  [  a  —  cos  (E  +  0)  )  ] 
xN  r  L  xN  J 


R  „  =  —  [  a  „  -  sin  (E  +  CD  )  ] 
yN  r  L  yN 
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Alternative  expressions  may  be  derived  for  use  in  special  circumstances, 
but  the  present  form  is  well  suited  to  the  low-eccentricity  orbits  of 
the  present  analysis. 

For  the  partial  derivatives  entering  into  r  Au,  the  differen¬ 
tials  of  rAv  are  obtained  from  equations  (8)  and  (9)  as  follows: 


Ar  cos  v  —  rAv  sin  v  =  Aa  (cos  E  —  e)  -  aAe  —  a  sin  EAE 


and 


Ar  sin  v  +  rAv  cos  v  =  Aa 


\fTT7  sin  E  -  *  e  stn  E  As 

yj  1  -  e2 


+  a  \l  1  —  e  cos  E  A  E 


from  which 


r  Av 


2 

e 


AE  + 


a  sin  E 


y 


1 


A  e 


Following  the  substitution  pattern  established  above  for  A  e  and  Ae 
leads  to 


Aa 


r  A v  =  V  AU  +V 

u  o  a  a 


+  V  M  Aa  M 
xN  xN 


v  M  Aa 
yN  yN 


(16) 


where  the  coefficients  or  partial  derivatives  are 


v  .d  /7T7 

u  r 

V  =  -  4  (U  -  U  )  V 
a  2  o  u 
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xN 


/ - 7  r  /7T7 

[  /  1  -  e  sin  E  cos  a)  (1  +  £)  +  —  ~e~ ~~ 


sin  CO  ] 


V  .  =  —  [  /  1  —  sin  E  sin  co  (1  +  — )  — 
yN  r  L  p'  e 


/l  -  e2 


COS 


0)] 


The  instability  of  these  derivatives  for  low-eccentricity  is  a  reflection 
of  the  indeterminate  nature  of  perigee  in  a  nearly  circular  orbit. 

The  rAoo  component  of  rAu  follows  directly  from  (13).  In 
combination  with  (16)  there  results 


rAu  =  r  Av  +  r  Ao) 


=  U 


AU  +  U 
o  a 


A  a 


-f  U 


xN 


A  a  M  +  U  M 
xN  yN 


Aa 


yN 


(17) 


with  the  coefficients 


u  =  V  =  —  \/  1  -  e2 
u  u  r 


U  =  V  = 
a  a 


-  -r  (U  -  U  )  U 
2  o'  u 


a2  .  „  r\  .  „  e2  -  (1  +  V  1  -  e2) 

U  x  =  —  !  (1  +  — )  sm  (E  +co)  +  a  „  e  sin  E  - * - L- 

xN  r  L  a  '  xN 


e  cos  E 


\/l  -  e2  (1  +  s/l  -  e2)2 


xN 


1  + 


n/T"-  e2 
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[  -  (1  +.  f)  cos  (E  +  oj)  +  ayN  e 


sin  E 


e“~  —  (1  -f  Vl  —  e'S  e  cos 

\Jl  -  e2  (1  +  yfl  -  e5) 


+ 


xN 


] 


1  + 


>4 


2 

e 


In  the  above  expressions ,  special  care  has  been  taken  in  the 
arrangement  of  terms  so  that  no  computational  problems  arise  as  the  eccen¬ 
tricity  approaches  zero.  These  expressions  for  Ar  and  r A  u  will  be 
combined  in  the  following  section  with  the  differential  expressions  in¬ 
volving  orbit  plane  orientation. 

b .  Differential  Expression  Extended  to  the  Orientation  of  the  Orbit  Plane 

The  orbit  plane  orientation  enters  into  the  analysis  through  the 
expression  for  the  observation  vector  p  ,  shown  in  Figure  III-2  and  defined 
by 


£-r+R-rU+R 

where  R  is  the  "station"  vector  from  observer  to  geocentric.  Differen¬ 
tiating  leads  to 


A£  -  AtD  +  r  AU  +  AR  .(18) 

For  the  present  analysis,  AR  is  assumed  zero,  and  no  attempt  to  correct 
station  location  will  be  made.  The  expression  for  A r  is  given  in  Eq. 
(15);  the  expression  for  AU  remains  to  be  determined. 

The  changes  in  U  can  be  attributed  to  two  sources;  the  change 
in  orientation  of  the  orbit  plane  and  the  change  in  position  of  the 
object  within  this  plane.  These  can  be  determined  from  the  identity 


U  =  N  cos  u  +  M  sin  u 


Au  =  A  N  cos  u  +  A  M  sin  u  +  Au  (M  cos  u  —  N  sin  u)  (19) 
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FIG.  III-2 .  OBSERVATIONAL  FRAMEWORK 
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The  quantities  of  A  N  and  AM  can  be  related  to  the  parameters  Aft  and  Ai 
through  the  following  line  of  reasoning.  Starting  with  the  definition  of 
the  orientation  parameters  in  terms  of  W 


W: 


W  =  sin  ft  sin  i 
x 


W  =  —  cos  ft  sin  i 

y 


W  =  cos  i 
z 


(20) 


Differentiation  leads  to 


AW  =  cos 

X 

ft 

sin 

i 

Aft 

+  sin 

ft 

cos  i 

A  i 

AW  =  sin 

ft 

sin 

i 

Aft 

—  cos 

ft 

cos  i 

A  i 

y 

AW  =  —  sin  i  A  i 
z 

The  coefficients  of  Aft  and  of  Ai  can  be  expressed  in  terms  of  M  and 
N  from  their  components. 


N: 


M: 


cos  ft 
sin  ft 
0 

—  sin  ft  cos  i 
cos  ft  cos  i 
sin  i 


(21) 


(22) 
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Thus 


AW  =  N  Aft  sin  i-M  Ai 


Similarly,  by  differentiating  (21)  and  (22)  expressions  for  A  N  and  AM 
follow. 

AN  =  M  Aft  cos  i-W  Aft  sin  i  (23) 

AM  =  —  N  Aft  cos  i  +  W  Ai  (24) 

By  combining  (18),  (19),  (23)  and  (24), the  expression  for  AjO  assumes 
the  form: 


A^  =  ArU  +  rA  u  (M  cos  u  —  N  sin  u) 


+  r Aft  [  cos  i  (M  cos  u  —  N  sin  u)  —  W  sin  i  cos  u)  ] 
4-  rA  i  (W  sin  u)  (25) 


or,  expressed  in  terms  of  the  U,  V,  W  unit  vectors 


A^  =  ArJJ  +  (rA  u  +  r  Aft  cos  i)  V 

+  (rAi  sin  u  —  r  Aft  sin  i  cos  u)  W  (26) 


The  six  orbit  parameters  enter  this  differential  expression  explicitly 
as  A  i  and  Aft  ,  and  implicitly  through  Ar  and  rA  u  as  given  by  (15) 
and  (17).  Then  combination  leads  to: 
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a£  -  AU0  111/*  uu)  *7(U1,U  Ua) 

+  a\n  <2rxn+2"xn>  +  a\n  «Vt!V  + 

+  A SI  (V  r  cos  i  —  W  r  sin  i  cos  u) 

+  A  i  (W  r  sin  u)  (27) 


where  the  partial  derivatives  denoted  by  R^,  Uu,  etc.,  have  been  pre¬ 
viously  defined. 

c •  Extension  to  Topocentric  Observation  Residuals 

The  final  step  in  the  derivation  of  the  differential  relation¬ 
ships  involves  the  translation  of  AjD,  as  expressed  in  (27),into  differ¬ 
ential  expressions  expressed  in  the  topocentric  observer1 s  geometry. 
Appropriate  observed  quantities  include  pairs  of  angles,  such  as  topo¬ 
centric  right  ascension  and  declination  or  altitude  and  azimuth,  and 
range.  These  observation  residuals  may  be  recognized  in  the  A £  expres¬ 
sion  by  writing 


A£  =  A(pL)  *  Ap  L  +  pA  L 


(28) 


and  the  definition  of  L  in  terms  of  topocentric  right  ascension  and 
declination: 


L: 


L  =  cos  6  cos  a 

X 


L  =  cos  6  sin  CX 

y 


L  =  sin  6 
z 


(29) 
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The  differential  AL  is  expressed  by 


AL:  - 


A  L  =  —  cos  j  sin  a  Aa  —  sin  5  cos  a  A6 

x 


AL  =  cos  5  cos  a  Aa  —  sin  6  sin  a  A6 

y 


AL  =  cos  6  A5 
z 


By  defining  two  additional  orthogonal  vectors,  A  and  D,  both  mutually 
perpendicular  to  L  and  with  A  lying  parallel  to  the  equatorial  plane 


A:  - 


D:  - 


A  =  —  sin  a 
x 


A  =  cos  a 

y 


A  =  0 
z 


D  =  —  sin  6  cos  a 
x 


D  =  —  sin  6  sin  a 

y 

D  =  cos  6 
z 


(30) 


(31) 


the  AL  can  be  expressed  as 


A  L  =  A  cos  6  Aa  +  D  A6 

or 

A£  =  L  Ap  +  A  p  cos  6  Aa  +  D  pA6  (32) 

Thus  the  successive  dot  products  of  A£  with  the  vector  triad  L,  A  and 
D  leads  to  scalar  differential  expressions  for  Ap  ,  p  cos  6  Aa  and  pA<5 
respectively.  In  practice,  the  vectors  L,  A  and  D  may  be  defined  either 
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by  the  observation  or  by  the  representation;  with  incomplete  observations, 
i.e.,  only  one  or  two  of  the  coordinates  such  as  obtained  by  range-only 
or  angle-only  instruments,  the  L,  A,  D  triad  must  be  obtained  from  the 
representation. 


Differential  expressions  involving  altitude  and  azimuth  obser¬ 
vations  njay  be  ^similarly  obtained  by  defining  the  horizon  system  vector 
triad  L,  A  and  D,  whose  horizon  components  are 


r 


—  cos  A  cos  h 

+  sin  A  cos  h 

+  sin  h 


(33) 


A: 


D: 


xh 


yh 


zh 


xh 


yh 


zh 


=  +  sin  A 

=  +  cos  A 

-  0 

■  +  cos  A  sin  h 

■  —  sin  A  sin  h 

■  +  cos  h 


(34) 


(35) 


where  A  is  ths  azimuth  angle  measured  east  from  north,  and  h  the  altitude 
angle  measured  from  the  local  horizon. 

These  horizon  components  must  be  rotated  into  the  equator  sys¬ 
tem,  to  which  the  components  of  A£  are  referred,  by 


Lx  "  Lxh  8x  +  Lyh  Ex  +  Lzh  Zx 


L  ■  L  .  S  +  L  .  E  +  L  .  Z 
y  xh  y  yh  y  zh  y 


Lz  -  Lxh  Sz  +  Lyh  Ez  +  Lzh  Zz 


L— »  A,  D 
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where 


S: 


E: 


Z: 


sin  0  cos  0 
sin  0  sin  0 

—  COS  0 

—  sin  0 
4-  cos  0 

0 

COS  0  cos  0 
cos  0  sin  0 
sin  0 


and  where  0  is  the  astronomical  latitude  and  9  is  the  local  sidereal 
t  ime . 

Then 


AL  =  A  cos  h  A  A  +  D  Ah 


and 


/\y 

A £  *=  L  A p  +  A  p  cos  h  A  A  +  D  Ah 


(36) 


Thus  scalar  differential  expressions  relating  the  orbit  parameters  to 
topocentric  altazimuth  observation  residuals  are  readily  obtained  by 
successive  dot  products  of  A^£  ,  given  by  (27),  with  L,  A  and  D  again 
determined  either  from  observation  or  representation. 
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APPENDIX  IV 


SCALAR  DIFFERENTIAL  EXPRESSIONS  FOR  DRAG 


The  development  of  the  secular  drag  perturbations  for  a  satellite 
moving  about  a  spherical  earth  has  received  considerable  attention  in  the 
literature  in  order  to  adapt  these  theories  to  a  differential- 

correction  procedure  for  earth  satellites  executing  grazing  entry,  where 
drag  represents  the  primary  perturbative  influence,  this  theory  has  been 
extended  to  include,  in  the  determination  of  density  variations  along  the 
path,  the  altitude  variations  arising  from  motion  over  an  oblate  earth. 

This  altitude  variation  is  significant  for  highly  inclined  satellites, 
amounting  to  21  km.  at  the  poles. 


IV.  1  DETERMINATION  OF  DENSITY  VARIATIONS  ALONG  PATH 

Altitude  variations  over  the  earth's  surface  arise  from  two  factors, 
as  illustrated  in  Fig.IV-1.  The  variation  due  to  eccentricity  is 

r  -  q  =  ae  (1  -  cos  E) 

from  the  perigee  distance  q.  A  second  altitude  variation  is  that  due  to 
oblateness,  where  the  earth  radius  varies  over  the  earth  as 

R  =  a  [  1  -  f  U2  +f  f 2  (U^-  U2)  ]  ,  U  *  sin  i  sin  (v  *4  gj) 
e  z  *  z  z  ’  z 

where  f  is  the  flattening  coefficient  (  =  1/298.3).  Combining  these  effects, 
the  altitude  variation  along  the  path  from  that  determined  at  the  two-body 
perigee  point  is  given,  to  first  order  in  flattening,  by 

ae  (1  -  cos  E)  +  f  (U2  -  P2  )  a 

z  e  e 
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FIGURE  IV- 1.  HEIGHT  ABOVE  OBLATE  EARTH 
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where  ^  is  determined  at  the  two-body  perigee  point.  This  expression  may 
be  expressed  in  terms  of  the  single  variable  E  by  the  substitution 

v  ■  E  +  e  sin  E  +  ... 

For  a  first  order  theory  in  flattening  and  eccentricity,  this  sequence  is 
terminated  at  v  *  E,  and 

U  *  sin  i  sin  (v+uu)  “  sin  i  sin  (E  +  (l)) 
z 

Further  manipulation  leads  to  an  altitude  variation  along  the  path  given  by 

ae  (1  -  cos  E)  +  —  8^n  ~  ag  [  cos  2u>  -  cos  2  (E  +  ou)  ] 

The  density  variations  along  the  path,  assuming  an  isothermal  (exponential) 
atmosphere  in  the  vicinity  of  perigee  altitude  is  then  given  by 

f  f  sin2i 

p  *  p  exp  <  -k  [  ae  (1  -  cos  E)  +  - - —  a  (cos  2cd  -  cos  2E  cos  2<jo 

"  1  2  *  ) 

+  sin  2E  sin  2cu  )  ]  > 

where  p^  is  the  perigee  density,  determined  from  an  appropriate  density  model, 
and  k~l  is  the  scale  height,  similarly  determined. 


IV. 2  SECULAR  VARIATION  IN  SEMI -MAJOR  AXIS 

Differentiation  of  the  vis-viva  (energy)  integral  leads  to  the 
expression  for  variations  in  semi-major  axis,  as  follows: 


2 

M* 

2  _ 

1  ' 

j± 

1  +  e  cos 

E  1 

S  “ 

r 

a  J 

a 

1  -  e  cos 

E 

and  2r  ■  £  -  2ss'  - 

~2 

or  s'* 

2a2 

.  A 

ss 

a 


where  r  is  the  velocity  vector  and  £  is  the  perturbative  acceleration. 


(1) 
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For  drag,  the  acceleration  component  s'  along  the  path  is  given  by 


1  Bp-  •  SSSS-* 

2  a  1  -  e  cos  E 


Following  the  procedure  outlined  in  ref.  (1),  the  variation  in  semi-major  axis, 
with  eccentric  anomaly,  is  derived: 

3/2  (2) 

(1  -  e  cos  E) 


as  given  in  ref.  (1).  Three  variable  expressions  are  involved,  and  the  following 
comments  are  germane: 

(T)  The  variations  in  density,  due  to  orbital  eccentricity  and/or 

oblateness,  are  very  large.  Eccentricity  effects  can  represent 
many  scale  heights  even  for  relatively  low  (-*0.01)  eccentricities, 
and  oblateness  can  lead  to  density  variations  of  the  order  of  one 
scale  height  for  highly  inclined  orbits. 

(T)  This  second  term  arises  from  variations  in  satellite  velocity,  and 
is  very  small  (~  3e)  for  low  eccentricity  orbits. 

(T)  This  third  term  arises  from  radial  distance  variations.  In  addition 
to  the  small  variation  (~  e)  for  orbits  of  Interest,  this  term 
partially  offsets  the  second  term  above. 


da 

dE 


r/a  v 

7T  a 


B  p  a 


6  (i/ 


[1  +  e  cos  E 

1  -  e  cos  E 


Thus,  the  significant  variations  in  semi-major  axis  arise  from  the  variations 
in  density  encountered  along  the  path.  In  the  following  analysis,  the 
secular  variation  in  semi-major  axis  shall  be  determined,  and  subsequently 
used  to  define  the  scalar  differential  expression  for  the  differential  correction 
of  drag. 

The  introduction  of  the  density  variations  along  the  path,  eqn.  (1), 
into  (2)  leads  to 


da  2 

—  *  -  Bp^  a  exp  (-kae)  <[  exp  (kae  cos  E)  •  (1  +  2e  cos  E  +..) 


exp 


2 

-kfa  *  (cos  2to  -cos  2u)  cos  2E  +  sin  2ou  sin  2E) 
e  2 
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The  last  exponential  term  is  expanded  in  a  Taylor  series  involving  terms  of  the 
form: 


a 

n! 


(1  -  cos  2E)  +  sin  2E  tan  2(JU 
with  the  observations: 


,  a 


kfsin^i 


a  cos  2(J u 
e 


(1)  Despite  the  factor  f  in  a,  a  is  not  a  small  number,  but  rather 
approaches  unity  at  altitudes  representative  of  terminal  decay. 

This  arises  from  the  factor  k,  the  reciprocal  scale  height,  and 
requires  that  sufficient  terms  be  carried  to  permit  the  factorial 
term  to  define  convergence. 

(2)  In  the  binomial  expansion  of  the  bracketed  quantity,  and  the  sub¬ 
sequent  integration  on  a  per-revolution  basis  to  define  the  secular 
variations  in  the  parameters,  all  odd  functions  of  E  (and  ou)  vanish. 


The  expansion  of  the  bracketed  quantity  leads  to  even  functions  of  E  of  the  form: 


cosrE  (1  -  cos  2E)8 

Their  subsequent  integration,  on  a  per-revolution  basis,  defines  a  sequence  of 
integrals  of  the  form: 

2tt 

f  exp  (z  cos  E)  F (cos  E,  1  -  cos  2E)  dE 
o 

where  z  =  kae 


These  integrals  define  Bessel  functions  of  imaginary  argument;  Table  IV-1  pro¬ 
vides  a  table  of  these  integrals  and  the  generating  series  for  the  B  (z)  terms, 
for  small  argument  (z  ^  2)  .  For  large  argument,  asymptotic  series  for  related 
functions  are  available. ^ 


By  performing  these  operations,  the  sequence  of  terms  for  the 
secular  variation  in  semi-major  axis  is  derived,  as  given  in  Table  IV-2  Terms 
of  order  e  have  not  been  included,  for  the  reasons  given  earlier. 
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TABLE  IV- 1 


BESSEL  FUNCTIONS  OF  IMAGINARY  ARGUMENT 
Integrals  of  the  Form 


2TT 


A  J 


2  tt 


exp(z  cos  E)  F(cos  E,  1-cos  2E)  dE 


Kernel, F 

Integral 

1 

B  (z) 

(l-cos-2E) 

o 

»1(z) 

(1-cos -2E)2 

3-B2(z) 

(1-cos  2E)3 

3-5-B3(z) 

(1-cos  2E)4 

3-5-7-B4(z) 

(1-cos  2E)5 

3-5-7-9-B5(z) 

2  cos  E 

z  B1(z) 

2cosE  (1-cos  2E) 

z  B2(z) 

2cosE  (1-cos  2E)2 

3-z  B3(z) 

2cosE  (1-cos  2E)3 

3 *5 *z  B4(z) 

2cosE  (1-cos  2E)4 

3-5-7-z  B5(z) 

2_ 
cos  E 

Bq(z)  -  h  Bl 

cos2E(1-cos  2E) 

B1(z)  -  h  [3-B. 

cos2E(1-cos  2E)2 

3-B2(z)  -  %  [3*5 

cos2e(1-cos  2E)3 

3-5’B3(z)  -  k  (3*5 

cos2E(1-cos  2E)4 

3-5-7  B4(z)  -  \  [3-5 

argument  (z  £  2): 

V‘>  - 

(z/2)Zr 

/  r!  (n+r)! 

r=0 
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TABLE  IV-2 


EXPRESSION  FOR  SECULAR  VARIATION  IN  SEMI -MAJOR  AXIS 


6a  2  ,  . 

5^  *  -Bpn  na  exp (-2) 


1 


2  3  4 

Bo  +  aBi  +  fp  (3>B2)  +  J\  <3'5’b3)  +  4T  (3"5'7*B4>  +  ••• 


+  (a  tan  2cd)' 


(B1~%  3-B2)  +  a(3*B2-%  3-5*B3)  +  jt  (3-5^-%  3'5’7-B^  +  .. 


+  -|  (a  tan  2ci))4 


(3-B2-3-5*B3  +  |  3'5*7-B4)  +  a(3  •  5^-3  *  5‘  7 +  £  3*5-7-9-B5)  + 


00 

o 


TABLE  IV-3 

EXPRESSION  FOR  SECULAR  VARIATION  IN  ECCENTRICITY 


eft  =  *Bpn  pn  exPC-z)fH' 


2  3  4 

B  +aB  +  ~  3‘B.  +  ^7  3'5'B  +  ^7  3'5'7'B  +  ... 

X  2  2  3  3 .  4  4  •  5 


+  (a  tan  2(d)' 


a 


(B2-^.3-B3)  +  a(3.B3-V3'5-B4)  +~  (3  •  S-B^-3 -5 •  7 -B^  + 


1  4 

+  —  (a  tan  2u)) 


(3-B3-  3-5*B4  +  £  3-5-7-B5)  +  a(3 *  5 *B4*3 -5 • 7 *B5  + 


•P'1 1- 


IV. 3  SECULAR  VARIATION  IN  ECCENTRICITY 

•k 

A  similar  procedure  leads  to  the  secular  variation  in  eccentricity, 
given  in  Table  IV. 3,  preceding  page. 

The  apparent  tangent  terms  in  Tables  IV  .2  and  IV  .3  present  no 
difficulty,  since  they  are  associated  with  factors  of  a  to  give: 

2 

a  tan  2cd  ■  -  kf  — ^ —  sin  2u) 

IV. 4  SCALAR  DIFFERENTIAL  EXPRESSIONS 

Secular  variations  in  semi-major  axis  and  eccentricity  have  been  pre¬ 
sented  above;  additional  secular  variations  in  inclination  perigee  and  node 
arise  from  atmospheric  rotation,  but  are  of  minor  importance  in  differential 
correction.  In  addition,  the  computational  investment  required  to  include 
short-period  terms  is  unwarranted  for  the  range  of  eccentricities  considered 
here. 


Orbit  parameters  included  in  the  differential  correction  include 
mean  motion  n,  a^  ■  e  cos  id  and  a  =  e  sin  u).  The  differential  correction 
for  drag  is  incorporated  by  expanding  the  scalar  differential  expressions 
(given  in  Appendix  I  for  the  two-body  parameters)  to  include  terms  in  AB, 

B 


as  follows,  where  Af  represents  an  observation  residual. 


An 


Af 


''An 

n 


3 

2 


1  6a 

—  *  T—  (t-t  ) 
a  o 


AB 

B  i 


+ 


.  a  I  6e 

+  xN - 

e  6t 


+ 


+ 


'Aa  „ 
y  N\ 


Aa  N  +  ayN  —  —  (t-t  ) 
y"  e  6t  ° 


CAO  AA  +  CAJ  Ai  +  CAll  AU 
An  Ai  AU  o 

o 


AB  | 
B 


See  Reference  (1)  for  partial  derivation. 
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6  3  6  0 

The  expressions  for  —  and  are  given  in  Tables  IV. 2  and  IV. 3.  The  apparent 

indeterminacy  of  the  and  coefficients  for  circular  orbits  is  more 

xN  yN 

6e 

apparent  than  real,  for  the  —  expression  contains  a  factor  z  =  kae.  In 

6 1 


practice,  the  e  term  in  z  is  associated  with  the  cos  u)  and  sin  cv  terms  to  avoid 
division  by  e. 
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APPENDIX  V 


SCALAR  DIFFERENTIAL  EXPRESSIONS  FOR 
SENSOR  CALIBRATION 


This  appendix  develops  the  4x7  matrix  of  coefficients  relating 
residuals  in  the  four  observation  components  -  range,  range-rate,  azimuth  and 
elevation  -  to  corrections  to  sensor  coordinates,  latitude,  longitude,  and 
height  above  the  reference  ellipsoid,  and  to  biases  in  range,  azimuth,  eleva¬ 
tion  and  station  time.  Although  the  determination  of  bias  in  range-rate  is  a 
trival  extension  of  this  theory,  the  existence  of  a  scale  factor  error  is  more 
probable,  and  is  revealed  in  a  characteristic  range-rate  residual  pattern. 

V.l  SCALAR  DIFFERENTIAL  EXPRESSIONS  FOR  SENSOR  LOCATION* 

The  problem  presented  is  that  of  correcting  a  geodetic  position  by 
using  observations  of  a  satellite  for  which  a  definitive  orbit  has  been  deter¬ 
mined.  To  elaborate,  assume  that  an  observation  network  has  been  established. 
The  stations  comprising  this  network  are  assumed  to  be  tied  together  exceedingly 
well;  that  is,  their  geodetic  positions  are  assumed  to  be  accurately  known. 
Further,  it  is  assumed  that  observations  from  the  stations  of  this  network  have 
resulted  in  the  determination  of  a  definitive  orbit  for  the  satellite.  Observa¬ 
tions  of  this  satellite  are  also  made  from  an  isolated  observing  site,  of  which 
the  approximate  location  is  known.  Corrections  to  the  site’s  assumed  position 


*  Koskela,  P.,  Nicola,  L.,  and  Walters,  L.  G.,  ’’Station  Coordinates  from 
Satellite  Observations,”  ARS  Journal,  Vol.  32,  p.  253,  February  1962. 
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must  now  be  obtained  by  using  the  satellite  observations.  The  position  obtained 
in  this  way  is  referred  directly  to  the  reference  spheroid  and  therefore  has  the 
advantage  of  being  independent  of  local  gravitational  anomalies. 

The  problem  may  be  visualized  with  the  aid  of  Figure  V-l.  As  the 
satellite  describes  its  orbit  about  Earth,  its  position  with  respect  to  the 
dynamical  center  is  defined  by  the  known  vector  r,  tabulated  as  a  function  of 
time.  The  observer's  position  is  given  by  R,  a  vector  conveniently  represented 
by  the  components  xc  and  y  in  the  meridian  plane  and  by  the  local  sidereal  time 
or  hour  angle  of  the  vernal  equinox  0.  For  a  fixed  observation  point  on  Earth, 
only  the  sidereal  time  varies  as  Earth  rotates. 

The  observations  are  components  of  the  vector £,  either  its  magnitude, 
slant  range  p ,  or  direction  cosines  L,  expressed  in  terms  of  the  topocentric 
angles  right  ascension  and  declination  or  elevation  angle  and  azimuth. 

The  station  vector  R  can  be  written  as 

R  *  -x  cos  0  I  -  x  sin  0  J  -  y  K  (1) 

—  c  —  c  —  c  ~ 


where  1 ,  J,  and  K  are  unit  vectors  directed  along  the  x,  y,  and  z  axes  respectively. 


9  ■  9gr  +  x: 


where  is  the  longitude  of  the  observer  measured  positively  to  the  east,  and 

0-,.,  is  the  Greenwich  sidereal  time.  0„  may  be  calculated  from 
GR  GR 

0GR  *  0GRq  +  0.004,375,26905  (t-tQ) 

where  0^  is  the  Greenwich  sidereal  time  in  radians  at  some  arbitrary  epoch, 
and  (t-t0?  is  the  mean  solar  time  interval,  in  minutes,  since  that  epoch.  The 
observer's  coordinates  in  the  meridian  plane  are 

x£  *  (C  +  H)  cos  <j) 
yc  *  (S  +  H)  sin 


where  (See  Section  V.3) 


North  Celestial  Pole 


Vernal 

Equinox 


FIGURE  V-l.  POSITION  RELATIONSHIP  OF  OBSERVER,  SATELLITE  AND  DYNAMICAL  CENTER 
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H  is  the  height  above  the  reference  ellipsoid  in  Earth  equatorial  radii,  and 


where  f  is  the  flattening  or  oblateness,  e  the  eccentricity  of  the  adopted 
reference  ellipsoid,  <j)  the  geodetic  latitude,  and  a^  Earth  fs  equatorial  radius. 

The  fundamental  relationship  between  the  positions  of  the  observer, 
dynamical  center,  and  satellite  may  now  be  written  as  (see  Figure  V-l) 

p  =  £  +  R 
or,  using  Eq.  1 

p  =  r  -  x  cos  9  I  -  x  sin  9  J  -  y  K 

This  equation  may  be  differentiated  to  relate  the  incremental  changes 
Ap,  residuals  in  range  or  angular  position  of  the  satellite  generally  taken  in 
the  sense  observed  minus  computed,  to  the  corresponding  changes  Ap  (orbit  un¬ 
certainties)  and  toAR  (station  errors).  If  it  is  assumed  that  the  orbit  is 
adequately  known,  it  will  remain  uncorrected  during  the  station  correction  pro¬ 
cedure,  and  the  incremental  changes  A  £  can  be  set  to  zero.  Following  this 
procedure 


Ap  =  Axq  [  -I_  cos  9  -  J  sin  9  ]  + 

AX  [  I  x  sin  9  -  J  x  cos  9  ] 

E  —  c  —  c 

This  differential  expression  may  be  transformed  into  the  more  familiar 
spherical  coordinates  involving  latitude  and  height  above  the  reference  ellipsoid; 
neglecting  second  and  higher  order  terms,  the  vector  differential  expression 
becomes 


Ap  =  [  yc  (!_  cos  9  +  J  sin  9)  -  x  K  ]  A^ 

+  [x  ( I  sin  9  -  J  cos  9)  ]  AX 

C  D 

A  H 

+  [ -  x  (I  cos  9  +  J  sin  9  )  -  y  K  ]  ---  „ 

L  c  —  —  Jc  —  C  +  H 
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An  alternative  expression,  in  terms  of  the  vector  triad  S  ,  E  and  Z  defining 
the  observers  south,  east  and  zenith  directions  may  be  written,  again  ignoring 
terms  of  order  e^  (flattening)  in  the  coefficients  for  this  differential 
expression: 

A£=  (C  +  H)  A<J)  S  -  (C  +  H)  cos  <|)AX.  e  -  Ah  z 

Successive  dot  products  of  this  vector  expression  with  the  triad  of  vectors 
L,  A^,  D^>  defined  by  the  observed  quantities,  leads  to  the  coefficients  for 
sensor  location  in  terms  of  range,  azimuth  and  elevation,  as  given  in  Table  V-l. 

The  comparable  coefficients  for  range-rate  may  be  derived  from 

pAp  =  =  £ •  A  R  ,  since  Ar  *  0. 

Taking  the  time  derivative, 

pA£  =  (£  -  £>  L)’A  R  +£•  A  R 

From 


R  =  -  (C  +  H)  Z 

AR  =  (C  +  H)  A<j>s  -  (C  +  H)  cos  X„A0E  -  AH  Z 


and  from 


These  substitutions  In  the  equation  forpA/b,  above,  lead  to 
pA/b  -  (C  +  H)  A<t>  [£*(£-  jbL  +  £  x  £)] 

-  (C  +  H)  cos<f)AXg  [  E*(£  -  pL  +  £x  Q)] 

-  Ah  [z  •(£_-  pL  +  £  x£)] 

These  expressions  are  given  in  Table  V-l. 

V.2  SCALAR  DIFFERENTIAL  EXPRESSIONS  FOR  TIME  AND  OBSERVATION  BIAS 

The  sensor  timing  error  At  may  be  identified  by  the  equivalent  orbital 
Au  (■  nAT  )  required  to  bring  the  orbit  timing  into  agreement  with  the  sensor 
observation°time8 .  These  coefficients  are  identical  to  those  employed  in  the 
orbit  differential  correction.  In  theory,  the  relation  of  timing  toAu  through 
the  mean  motion  n  is  strictly  valid  for  low-eccentricity  orbits;  in  practice, 
the  iterative  application  of  these  expressions  produce  rapid  convergence. 

The  additional  terms  in  Table  V-l relating  to  biases  in  the  observed 
quantities,  are  trivial. 

V.3  COORDINATES  OF  THE  SENSOR  AND  THE  SATELLITE 


This  section  is  appended  to  derive  the  expressions  for  sensor  location 
used  in  the  program  and  the  expression  for  the  height  of  a  satellite  above  the 
oblate  spheroid.  The  last  expression  is  cited  on  Page  4  of  Appendix  I  and 
Page  1  of  Appendix  IV. 


The  meridian  section  of  the  earth  is  assumed  to  be  an  ellipse.  The 
coordinates  of  the  intersection  of  the  radius  vector  with  the  spheroid  are 


i 


cos 

0  = 

a  cos  E 
e 

sin 

1! 

a  Vl-e2 

sin  E 
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TABLE  V-l.  SCALAR  DIFFERENTIAL  EXPRESSIONS  FOR  SENSOR  CALIBRATION 


-  SENSOR  LOCATION  ERROR  - 

(C+H)  A <f>  -(C+H)  cos  0 AX  -AH 


TIMING 
n  AT 


-  OBSERVATION  BIASES- 


B 


vo 

o 


RANGE 

Ap 


AZIMUTH 
p  cos  h  A  A 


ELEVATION 

pAh 


S-4 

1 

0 

0 

Z-^  (-0) 

4 '  ©>u  +  ™u> 

0 

P  cos  h 

0 

E 

- n 

Z-Dh 

D. • (UR  +  VU  ) 

— h  —  u  —  u 

0 

0 

p 

S-(^  -  /5>L  + 

E-(£  -  /&L  +  £^1) 

Z  (£_  -  Pl  +  jQxft) 

(A  -  PL) * (URu  +  VUu) 

+  £  [U(R  +  vU  )  + 

—  u  u 

V (0  +  vR  )] 

—  u  u 

0 

0 

0 

RANGE -RATE 
pAf) 


where  x^9  y^,  a^  and  e  are  the  same  as  in  V.l  above  and  0fis  the  geocentric 

latitude,  E,  the  eccentric  (or  reduced)  latitude,  corresponds  to  the 
eccentric  anomaly  in  the  orbital  ellipse. 

2.  2  ,  2  2,-2  2  ,-iv  2  ,-  2N 

(1-e  )  x  +  y  =  r  (1-e  cos  0  )  =  a  (1-e  ) 
c  c  c  e 

a.  Geodetic  Sensor  Location 


Consider  a  small  change  of  eccentric  latitude.  This  will  produce 
(see  Figure  V-2)  the  coordinate  changes 

-dx  =  ds  sin  0  =  a  dE  sin  E 
c  e 


~\j  2 

dy  =  ds  cos  0  =  a  dE  *l-e  cos  E 
J  c  e 


where  ds  is  the  displacement  along  the  meridian. 
Thus, 


(1-e2)  (dxc)2+(dyc)2=(ds)2(l-e2sin20)  =  (aedE)2(l-e2) 


from  which 


.  _  ds  ,  -  _  V 1-e"  sin  0 

sin  E  =  - — —  sin  0  =  - e 


Vl-e* 


a  dE 
e 


Vl-e2  sin^  0 


and 


V 1-e2 cos  E  =  da 


A  Vl-e2  cos  0 

nr  cos  0 

e  Yl-e  sin  0 


For  a  point  on  the  spheroid  then 

a  cos  0 


=  C  cos  0 


and 


a  (1-e2)  sin  0 
e 

Vl-e^  sin2  0 


=  S  sin 
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ro 
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FIGURE  V-2.  MERIDIAN  CROSS-SECTION  OF  EARTH 


Since  height,  H,  Is  ordinarily  measured  opposite  to  the  direction  of  gravity, 
for  a  general  point  and  In  particular  for  the  sensor 

x  *  (C  +  H)  cos  0 
c 

yc  =  (S  +  H)  sin  0 

as  in  V.l  above. 

b .  Height  Above  Oblate  Spheroid 

In  terms  of  geodetic  latitude,  the  distance  from  the  center  of  the 

earth  is 


r 


c 


r 


c 


Vi  2~  2\  ,  2  .  2  z 

=  * 1-e  (1  -  e  )  +  e  sin  0 


r 


c 


a 


e 


since 


U  =  sin  0 
z 
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APPENDIX  VI 


VARIATION- OF -PARAMETERS  FOR 
NON  EQUATORIAL  ORBITS 


In  the  variation-of -parameters  technique,  the  parameters  of 
an  osculating  orbit  are  determined  by  the  actual  position  and  velocity 
at  a  given  instant.  These  parameters  describe  the  two-body  orbit  that 
the  object  would  follow  if  all  subsequent  perturbations  were  removed. 
Thus,  in  Fig.VI-l,the  full  curve  represents  the  actual  disturbed  path 
of  an  object,  and  the  two  dotted  ellipses  represent  the  osculating  two- 
body  orbits  derived  from  position  and  velocity  at  the  points  A  and  B, 
respectively . 

When  the  two  points  are  close  together,  the  elements  of  the 
two  osculating  ellipses  will  differ  by  relatively  small  amounts. 

Hence,  these  elements  may  be  visualized  as  parameters  that  vary  with 
time  because  of  the  disturbing  forces. 

The  ephemeris  program  employs  the  variation-of -parameters 
technique  and  thereby  overcomes  the  computational  inefficiency  of  nu¬ 
merical  Integration  methods  where  the  total  acceleration  is  computed 
and  integrated  step-by-step.  Computational  efficiency  is  achieved  by 
replacing  the  integrated  variables  by  parameters  which,  in  the  two-body 
case,  are  invariant.  Thus,  the  numerical  integration  required  relates 
only  to  the  components  of  the  force  field  which  disturb  the  two-body 
motion.  These  perturbations  are  extremely  small  in  comparison  to  the 
central  force  term  during  normal  orbiting  flight.  An  alternative  per¬ 
turbation  technique,  the  so  called  "Encke"  method  which  integrates  de¬ 
partures  from  the  two-body  orbit,  is  not  well  adapted  to  the  integra¬ 
tion  of  a  satellite  orbit  when  the  object  is  exposed  to  continuous  drag 
over  the  orbit. 
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In  the  following  analysis,  the  time  rate -of -change  or  varia¬ 
tion  of  functions  will  be  divided  into  two  parts  as  follows: 

If  f  is  any  function, 

i_4£.d£sf+fN  (1) 

ke  dt  dT 


where  T  =  k0  (t-tQ),  the  normalized  time. 

In  Eq.  (1)  f  (f-dot)  is  the  "two-body  variation"  that  would 
hold  at  the  instant  if  all  subsequent  perturbations  were  suddenly  re¬ 
moved,  i.e.,  the  variation  in  the  instantaneous  osculating  orbit;  and 
f  (f -grave)  is  the  perturbative  variation,  i.e.,  the  part  of  the  varia 
tion  caused  by  the  disturbing  forces  (in  this  case  the  atmospheric  drag 
and  Earth's  bulge). 

For  some  functions  the  perturbative  variations  are  zero  be¬ 
cause  the  instantaneous  osculating  orbit  is  defined  so  as  to  yield  for 
their  variations  the  same  values  as  the  actual  orbit.  Thus,  for  x,  y, 
and  z  referred  to  non-rotating  axes, 


dx 
d  T 


=  0 


x 


* 


y>z 


(2) 


and  for  velocity  functions  independent  of  orientation, 


dr 

dr 


r'  =  0 


(3) 


♦Denotes  that  similar  relations  hold  for  y  and  z  components. 
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ds  = 

dr 


s' 


=  0 


(4) 


For  the  elements  that  would  be  constant  in  two-body  motion, 
the  variations  are  consequences  of  the  disturbing  forces  alone.  For 
example , 


da 

dr 

dn 

dr 


a' 


n'  , 


a  =  0 


n  =  0 


(5) 

(6) 


The  anomalies,  v  and  E,  and  the  mean  longitude,  L,  have  both 
kinds  of  variation  because  they  are  referred  to  a  perturbed  or  "accel- 


origin. 

Thus , 

dv  _ 

dr 

v  +  vs  , 

If 

CM 

1 

U 

II 

•  > 

( 7a, b) 

di  = 
dr 

E  +  E'  , 

£  =  r“*  s/  Jl  /a 

(8a, b) 

dL  x 
dr 

i  +  i 

(  9a) 

Since  L  =  M  =  n/ke,  the  last  of  these  expressions  is  more  familiar  in 
the  form 


—  =  n  +  k  t 

dt  e 


(9b) 


The  parameters  which  are  integrated  in  the  low-eccentricity 
ephemeris  program  are: 


h  =  \Tp  W 
a  =  e  P 


the  orbital  angular  momentum  vector 
per  unit  mass 

-  a  vector  defining  eccentricity  and 
perigee  location 
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L  =  M  +  0)  +  -  the  mean  longitude  of  the  object 


The  vector,  h,  is  a  logical  choice  for  satellite  orbits 
since  it  embodies  the  orientation  of  the  orbit  plane  as  well  as  the 
two-body  properties  of  the  orbit.  The  vector,  a^  combines  the  effects 
of  eccentricity  and  perigee  location  in  a  way  which  maintains  continui¬ 
ty  through  zero  eccentricity.  The  remaining  parameters  define  the  po¬ 
sition  of  the  object  in  its  orbit.  Seven  quantities  are  integrated, 
including  six  components  of  the  two  vectors,  rather  than  the  minimum 
six,  to  avoid  complications  in  the  evaluation  of 

some  of  the  parameters. 


In  the  application  of  the  variati  on-of -parameters  technique, 
the  perturbative  variations  in  the  parameters  must  be  related  to  the 
components  of  the  disturbing  accelerations.  The  perturbative  variation 
in  n,  the  mean  motion,  follows  from  the  definition, 


n 


leading  to 


rt 


ke  ^ 


3 

2 


n  a  ^ 


(10) 


(ID 


To  relate  n  to  the  force  field,  it  is  necessary  to  express  &  in  terms 
of  the  perturbative  variations  in  velocity,  This  is  readily  accom¬ 

plished  through  the  vis-viva  integral, 


.2 

s 


r  -r 


=  Mt  —  ■  —  ] 


r  a 


(12) 
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Differentiation  yields 


2r 


dr 
d  t 


24  dr  (  M  da 
r2  dT  a 2  dT 


(13) 


However,  r  is  independent  of  orientation,  and  the  semi-major  axis  and 
mean  motion  are  constants  in  two-body  motion.  Therefore,  using  Eqs. 

(3)  and  (5),  and  the  definition  given  in  Eq.  (1)»  leads  to 

2r  •  (r  +  £  )  =  -  2/£  r  +  iL  a  (14) 

rz  a2 

The  first  terms  on  each  side  of  Eq.  (14)  are  equivalent,  since  they 
arise  in  the  differentiation  of  the  two-body  problem.  The  perturbative 
component  of  a  is  therefore 

i  =  2air-r  (15) 

Combining  Eq.  (15)  with  Eq.  (11)  yields 

ju  li  =  -  3na  r/£  Q.6) 

Eq.  (16)  then  relates  i)  to  the  forces  tending  to  disturb  the  two- 
body  motion.  These  forces  are  embodied  in  £  and  include  contributions 
from  drag  and  the  Earth's  bulge.  It  is  interesting  to  note  that  the 
mean  motion  is  affected  only  by  the  tangential  component  of  r  ;  drag 
is,  therefore,  responsible  for  the  secular  variations  in  n. 

By  similar  reasoning,  other  parameters  may  be  related  to  the 
non  two-body  components  of  the  force  field.  The  angular  momentum,  h, 
is  immediately  expressed  through  its  definition, 

*//7  h  =  r  x  r  (  17) 


Therefore,  by  a  shortcut  differentiation  process  that  may  be  inferred 
from  the  foregoing  differentiation  of  r/£, 


99 


(18) 


JJi  ti  =  r  x  £ 


because  h  and  i*  are  both  zero. 

Next,  the  vector,  a,  can  be  expressed  as 

*/jU  a  =  6  r  -  D  r 

where  D  is  defined  as  follows: 

•Jj!  D  =  r  •£  and  s/ji x  D  =  r  •  r_  - 

r 

Differentiation  leads  to  the  expression  for  , 


(19) 


(20) 


m£  =  (2L-l  )  l  ■  <l-£  )  i  -  (i-£)  £ 


(21) 


The  perturbative  expression  for  the  mean  longitude,  L,  is 
complicated  by  its  dependence  upon  the  rotation  of  the  orbit  plane 
about  the  instantaneous  radius  vector,  r,  as  well  as  the  transverse  ac¬ 
celeration  components.  Thus, 


i)  =  Mv  +  (22) 

for  direct  motion. 


While  each  of  the  three  angles  on  the  right  - 
nand  side  in  Eq .  (22)  behaves  badly  for  small  eccentricity  or 

small  inclination,  their  sum  ii  is  well  behaved. 

The  vectors,  r_  and  U,  are  not  instantaneously  affected  by  the 
perturbations,  so  Harises  solely  from  the  perturbations  of  the  orbital 
axes.  From  Fig.  VI-2(a) 

^  =  p-{£  =  -  &■£ 
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but  a  =  e  P 

so  -evx  =  =  Qxax  Qy^y  Qz^z  (23) 

The  orthogonal  component  of  the  perturbative  acceleration  is 

r£  =  W-r  =  Wxx  +  Wyy  +  Wzz  (24) 

where  W  is  the  unit  vector  perpendicular  to  the  instantaneous  orbit 
plane . 

The  perturbations  of  the  orientation  angles  are  caused  by  the 
rotation  of  the  orbit  plane  about  the  instantaneous  radius  vector,  r. 
From  Fig.  VI-2(c) 

N  =  U  cos  u  -  V  sin  u  ,  (25) 

and 

M  =  U  sinu  +  Vcos  u,  (26) 

where 

u  =  v  +  a)  «  (27) 

From  Fig.  VI-2(b),  looking  along  the  instantaneous  line  of  nodes, 

uN  =  -  rf  cos  i  (28) 

Now,  at  any  instant,  the  true  longitude  is 

H  *  u  +  Cl  ( 29) 

ao  that 

/  ■  u'  +  ft  -  rt  (1-cot  i)  (30) 
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and  the 


But  from  Fig. VI-2 (b),  Eq.(25), 

U-W'=  -W*U'=  0, 

sin  i  =  -W-tf  =  N-tf 


fact  that  =  0,  so  that 

-V  V  sin  u  (31) 


Now,  V  =  WxUso’VpV  =  hxU  =  r  x  r  x  U 
or  ^pV  =  (rr  -  rr)  /  >Jjj. 

so  that 

\Tp  V  +  =  (rr  -  rr)  /  4\±  (32) 

because  rN  =  0  and  =  0 . 

Now  dotting  W  into  Eq .  02)  above  and  noting  that 

W-V  =  0  , 
v/*v  =  -w-v' 

and 

W-r  =  0 

because  of  orthogonality, 

V/W  =  r  =  -V-w'  ,  (33) 

Mjd  p 

Substituting  Eq .  (33)  into  Eq.  (31)  and  recalling  Eq.  (24), 

St  =  L2j-Si.n  u  .  (34) 

sin  i 

Finally,  entering  these  relationships  into  Eq.  (30), 

=  r  s^-n  u  (l~cos  i)  rl? 

•T\i  p  sin  i 
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I 


or 


/  =  (r  sin  u  sin  1)  rb 
J~iT p  (1+cos  i) 


But 


and 


•o 


z  ■  r  sin  u  sin  i 

cos  i  ■  , 

a  z  (rb) 

t  ■  -  . 

(14«  ) 


(35) 


ing  t  : 


The  following  fundamental  relationships  are  used  in  dariv- 


M  -  E  -  e  sin  E  ,  (36) 

e  sin  E  -  -£U  9  (37) 

r  cos  v  *  a  (cos  E  -  e)  ,  (38) 

r  sin  v  *  a  -eZ  sin  E  ,  (39) 

j  -•  (1-e  cos  E),  (40) 

l  -  M+oo  +  n-  M+  n,  (4i) 

n  -  <i)  +  n  .  (42) 


Differentiating  Eq.  &2) , 

If  -  0)  +  ft  -  £S  -  v'. 


(43) 


so  that 


I)  -  M'+  fr  -  ^ 


(44) 
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[ 


Differentiating  Eq .  (40)  and  Eq .  (37), 


(e  cos  e)  =  eN  cos  E  -  e  EN  sin  E  =  — 9  av  , 


(45) 


1  aN 


(e  sin  E^  =  ex  sin  E  +  e  EN  cos  E  =  ,■  —  —  —  e  sin  E,  (46) 

n/  a  2  a 


from  which 


rr 


jin  E  +  —  —  f  (1+  1)  cos  E  ■  e] 

VTTi  2  a  1  a  ’ 


(47) 


and 


e  Ev  =  —  cos  E  -  i-  —  (§■  +  1)  sin  E 


(48) 


Now  differentiating  Eq.  (38)  and  rearranging  terras, 


a  e'  sin  E  =  rv'  sin  v  +  aN  (cos  E  -  e)  -  ae'  . 


(49) 


Substituting  Eq .  (39)  and  Eq .  (40)  and  dividing  out  a  sin  E  , 


E'  =  v'VT^2  -  _e£,  -  1  i.  e  sin  E  . 

N/7Ti  2  a 


(50) 


Differentiating  Eq .  (36)  and  introducing  Eq .  (46), 


M  =  E  -  (e  sin  EJ  =  E‘  - 


rr 


+  —  —  e  sin  E 
jj.  a  2  a 


(51) 


Combining  Eq.  (50)  and  Eq .  (51), 


Mv  =  ^  ^7e2  -  -jUi  . 

vju  a 

Finally,  using  Eq .  (44)  and  Eq.  (52), 


(52) 


l)  =  £'  - 


2  rr 


■s/juTa"  1  +  *J  1-e2 


(53) 
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p-(£ 


-  2 :t 
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APPENDIX  VII 


LIST  OF  SYMBOLS 


Many  of  the  symbols  used  in  this  document  occur  with  the  same 
meaning  throughout  the  various  volumes.  These  have  been  defined  below. 
We  have  attempted  to  avoid  inventing  new  notation  wherever  possible  . 

The  symbology  is  closest  to  that  of  Baker  and  Makemson,  Introduction 
to  Astrodynamics ,  Academic  Press,  1960.  This  in  turn  is  based  on  the 
notation  of  Herrick, who  follows  closely  the  recommendations  of  the 
International  Astronomical  Union. 

Some  symbols, which  are  used  only  once,  are  defined  at  that 
place  and  not  in  this  list.  Subscripts  are  sometimes  an  obvious  use 
of  a  symbol  defined  in  this  list  to  delimit  another  symbol.  Such  a 
subscript  is  not  generally  defined  separately.  For  any  other  omissions 
apologies  are  offered. 

Vectors  are  denoted  by  underlines.  While  most  vectors  symbol¬ 
ized  with  capital  letters  are  unit  vectors,  R,  for  instance,  is  not. 


Superscripts 
•  (Dot) 

1  (Grave) 

-  (Bar) 


Derivative  with  respect  to  T  (canonical  time) 
Perturbative  derivative  with  respect  to 
Average 
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LIST  OF  SYMBOLS 


Azimuth  angle,  measured  eastward  from  north  in  the  horizon 
plane 

Effective  cross-sectional  area  of  the  satellite 

Unit  vector  perpendicular  to  the  line-of-sight  vector, 

L,  in  the  direction  of  increasing  right  ascension 

Unit  vector  perpendicular  to  the  line-of-sight  vector^L, 
in  the  direction  of  increasing  azimuth 

Geomagnetic  index 

Semi-axis  major  of  ellipse  or  hyperbola,  mean  distance 
Mean  equatorial  radius  of  the  earth 
eP,  vector  integral  of  the  equations  of  motion 
e  cos  x  =  a  •  N 
e  sin  jo  =  a  •  M 

Subscript  denoting  bulge  (asphericity)  perturbations 

cda 

Ballistic  parameter,  — jjj- 

Modified  Bessel  function  of  imaginery  argument,  Z: 


2  n 


I  (Z) 
n 


z 


Bias  in  0.  stored  in  BIBUF 


l 


Semi-axis  minor:  b 


a 


Auxiliary  quantity  used  to  find  station  coordinates: 


C 


sin^  0 


Total  atmospheric  drag  coefficient 

Coefficient  of  ^  i  in  the  equations  of  condition  (differential 
correction),  where  i  is  a  parameter  to  which  correction  is 
sought 

Coefficient  of  cos  j  X  in  tesseral  harmonic  of  order  i 


Speed  of  light 

Subscript  indicating  computed  value  or  a  quantity  derived 
from  the  ephemeris  computation 


Auxiliary  variable:  ^ (X  D  =  rr 
Day  Number 

Subscript  denoting  drag  (atmospheric)  perturbations 

Unit  vector  perpendicular  to  line-of-sight  vector,  L, 
in  direction  of  increasing  declination 

Unit  vector  perpendicular  to  line-of-sight  vectors,  L, 
in  direction  of  increasing  elevation  angle 


Differential  operator 


Eccentric  anomaly 

Unit  vector  to  east  point  of  horizon 


Flux  of  solar  radiation  at  10.7  cm 


F10  averaged  over  three  solar  rotations 
Magnitude  of  the  force  of  solar  radiation  pressure 

Flattening  of  earth  spheroid 

Acceleration  of  gravity 

Components  of  the  perturbative  acceleration  due  to  the 
earth's  asphericity  in  the  easterly,  southerly  and  radial 
directions.  Note  that  radial  is  defined  here  in  the 
geocentric,  not  geodetic,  sense  and  thus  southerly  is 
perpendicular  to  the  geocentric  radius  and  not  in  the 
local  horizon  plane 

Subscript  denoting  Greenwich  meridian 

Scale  Height 

Unit  vector  perpendicular  to  line  of  nodes  in  the  equator 
plane 

Hour  angle  of  the  sun 
Density  scale  height 

Height  above  sea  level 

Subscript  indicating  horizon  coordinate  system 
Elevation  angle 

Vector  integral  of  the  equations  of  motion:  yr  w 


h  is  the  angular  momentum  per  unit  mass 
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In(z) 


Unit  vector  toward  the  foot  of  the  mean  vernal  equinoctial 
colure  of  epoch  on  the  true  equator  of  date 

Bessel  Function  of  imaginary  argument  and  order  n 


inclination  of  the  satellite  orbit  to  the  equator  plane 


Unit  vector  to  foot  of  summer  solstitial  colure  on  the 
true  equator  of  date 


K 


Unit  vector  to  the  north  pole  of  the  true  equator  of  date 


k 

k 


reciprocal  of  scale  height 

Gaussian  gravitational  constant,  square  root  of  the  product 
of  the  universal  constant  of  gravitation  and  the  mass  of 
the  earth,  k^  =  0.07436662  (earth  radii)3/2  per  minute 

Dimensionless  number  equal  to  ke,  ratio  of  canonical  to 
ordinary  solar  time. 


L 

L 


Mean  orbital  longitude 
Line-of-sight  unit  vector 


True  orbital  longitude 


M 

M 


Mean  anomaly 

Unit  Vector  toward  northern  antinode  of  the  orbit 
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m 


mass  of  satellite 


N  Unit  vector  toward  ascending  node 


n  Mean  angular  motion  of  satellite 


0^  General  notation  for  an  observed  quantity 

o  Subscript  referring  to  epoch,  tQ 

o  Subscript  indicating  observed  value  or  a  quantity  derived 

from  observed  values 

0  Subscript  denoting  sun 


P 

P 


P 


P 

n 


P 

n,m 


P' 

n 


P' 

n,m 


pj 


Period  of  revolution 
Pressure 

Unit  vector  to  perigee 
Solar  radiation  pressure 
Legendre  polynomial  order  n 
Associated  Legendre  function 

Derivative  of  the  Legendre  polynomial  with  respect  to 
its  argument 

Derivative  of  the  associated  Legendre  function  with 
respect  to  its  argument 

Semi-latus  rectum,  parameter  of  conic  section 

General  notation  for  a  parameter  describing  the  satellite 
motion 
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Q  apogee  distance 

£  Unit  vector  parallel  to  velocity  at  perigee 

q  perifocal  distance 


R 

R 

R 


n,m 


Distance  of  observing  station  from  center  of  earth 

Subscript  denoting  radiation  pressure  perturbations 

Vector  from  observing  station  to  center  of  earth.  Its 
components  are  X,  Y ,  Z 


P 

n,m 


d  -  uz 


2^—1/ 2 


R'  P'  (1  -  U  2)1/2 

n,m  n,m  z 


Radial  distance  from  center  of  earth  to  satellite 
Position  vector  of  satellite.  Its  components  are  x,  y,  z 
Correlation  between  corrections  to  elements  i  and  j 


S 


2 

C(1  -  e)  ,  auxiliary  quantity  used  in  finding  station 
coordinates 


Unit  Vector  to  South  point  of  horizon  or  similar  unit 
vector  perpendicular  to  radius  vector. 

Coefficient  of  sin  j  X  in  tesseral  harmonic  of  order  i 
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Fraction  of  interpolation  interval 


arc  length 

speed,  magnitude  of  velocity  vector 


Temperature 


time 


Mean  argument  of  latitude 
Unit  vector  toward  satellite 


True  argument  of  latitude 


Unit  vector  in  orbit  plane  perpendicular  to  radius  vector 
in  direction  of  increasing  anomaly,  transverse  unit  vector 


true  anomaly 


Unit  vector  perpendicular  to  orbit  plane  in  same  direction 
as  angular  momentum  vector 


R  •  1^,  component  of  station  vector 


£  •  I ,  component  of  satellite  position  vector 


J  •  R  component  of  station  vector 


J  •  i:  component  of  satellite  position  vector 


K  •  R  component  of  station  vector 


K  •  r^  component  of  satellite  position  vector 
ae 

~  argument  in  the  analytic  integration  of  drag 
P 


right  ascension 
2 

l/2kf  sin  i  a  cos  2  w  5  argument  in  the  expansion  of  the 

analytical  Hrag  perturbations  to  account  for  the  oblateness 
of  the  atmosphere 


Reflectivity  of  satellite 


Differential  operator  producing  a  small  finite  increment, 
especially  a  residual 


Declination 


Small  criterion  for  convergence  or  error  test 


Obliquity  of  the  ecliptic 


sin  *  (r  *),  auxiliary  angle  used  in  radiation  pressure 
perturbations 

Minimum  angular  distance  from  subsolar  bulge  (Appendix  II, 
Volume  I) 


Sidereal  time 


Upper  bound  (fraction)  for  correction  to  ballistic  parameter 

Tabular  coefficient  of  ballistic  parameter  to  account  for 
variation  of  drag  coefficient  with  altitude 


Geographic  longitude  measured  east  from  Greenwich 
Longitude  west  of  Greenwich 

Gravitation  constant,  usually  unity 

speed  with  respect  to  atmosphere,  ground  speed 
Velocity  with  respect  to  atmosphere,  ground  velocity 

Maximum  angular  distance  from  subsolar  bulge  (Appendix  II, 
Volume  I) 

Longitude  of  perigee 
subscript  denoting  perigee  value 


Atmospheric  density 

Range,  distance  from  sensor  to  satellite 
Position  vector  of  satellite  with  respect  to  sensor 

Summation  operator 

Standard  deviation 

canonical  time  r  =  k'(t  -  t  ) 

o 

Modified  hour  angle  of  subsolar  bulge  (Appendix  II, 
Volume  I) 

Potential 

Geodetic  or  astronomical  latitude 
Geocentric  latitude 

Elongation  from  Sun 

Right  ascension  of  the  ascending  node 
Argument  of  perigee 
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